Changes in the serum levels of IL-1 family cytokines after kidney allograft transplantation
Statistical report - data exploration and statistical analysis
Authors’ affiliations:
Cecrdlova E1, Krupickova L1, Fialova M1, Novotny M2, Svachova V1, Mezerova K1, Viklicky O2, Tichanek 3, Striz I1
1Department of Clinical and Transplant Immunology
2Department of Nephrology
3Department of Data Science
Institute for Clinical and Experimental Medicine, Prague, Czech Republic
The paper is under review at the Cytokine journal
1 Packages and functions
Open code
if (T) {rm(list = ls() )}
if (T) {
suppressWarnings(suppressMessages({
library(tidyverse)
library(stringr)
library(ggpubr)
library(emmeans)
library(gtsummary)
library(car)
library(sjPlot)
library(flextable)
library(openxlsx)
library(mgcv)
library(pROC)
library(cowplot)
library(boot)
library(glmnet)
library(brms)
library(projpred)
library(RJDBC)
library(janitor)
library(arm)
library(corrplot)
library(lubridate)
library(kableExtra)
library(ggdist)
library(bayesplot)
library(coxed)
library(quantreg)
library(rstudioapi)
library(effects)
library(lqmm)
# Functions clashes
select <- dplyr::select
rename <- dplyr::rename
mutate <- dplyr::mutate
recode <- dplyr::recode
summarize <- dplyr::summarize
count <- dplyr::count
# Simple math functions
logit <-function(x){log(x/(1-x))}
inv.logit <- function(x){exp(x)/(1+exp(x))}
}))
}1.1 run_model function
The function to (i) load OR (ii) run & save (if has not been done previously) results of any computation or complex table production
Open code
run <- function(expr, path, reuse = TRUE) {
# Initialize 'fit' to NULL to ensure it's always defined
fit <- NULL
# Construct the file path only if 'reuse' is TRUE
if (reuse) {
path <- paste0(path, ".Rds")
fit <- suppressWarnings(try(readRDS(path), silent = TRUE))
# Check if 'fit' is an error (file not found or couldn't be read)
if (inherits(fit, "try-error")) {
fit <- NULL
}
}
# If 'fit' is NULL (either because 'reuse' is FALSE, or the file couldn't be read), evaluate 'expr'
if (is.null(fit)) {
fit <- eval(substitute(expr))
# Save 'fit' only if 'reuse' is TRUE and a valid 'path' is provided
if (reuse && !is.null(path) && nzchar(path)) {
saveRDS(fit, file = path)
}
}
return(fit)
}1.2 Data upload
Open code
patient_table <- read.table('data/pat_data.txt')
patient_table <- patient_table %>% select(-c(Patient, date_of_rejection,
bid, Dg))
summary(patient_table)
## id rejection_afterTX group receiver_sex
## Length:186 Min. : 6.00 Length:186 Length:186
## Class :character 1st Qu.: 17.25 Class :character Class :character
## Mode :character Median : 69.50 Mode :character Mode :character
## Mean :100.30
## 3rd Qu.:110.25
## Max. :505.00
## NA's :142
## receiver_age donor_sex donor_age N_mismatch
## Min. :22.00 Length:186 Min. : 1.00 Min. :0.000
## 1st Qu.:45.00 Class :character 1st Qu.:44.00 1st Qu.:2.250
## Median :56.50 Mode :character Median :55.00 Median :3.000
## Mean :54.67 Mean :52.80 Mean :3.304
## 3rd Qu.:65.00 3rd Qu.:64.75 3rd Qu.:4.000
## Max. :80.00 Max. :81.00 Max. :6.000
## NA's :48 NA's :48 NA's :48
## dsa anti_hla PRA_max PRA_actual
## Length:186 Length:186 Min. : 0.00 Min. : 0.0
## Class :character Class :character 1st Qu.: 2.25 1st Qu.: 0.0
## Mode :character Mode :character Median :10.00 Median : 2.0
## Mean :19.09 Mean :11.6
## 3rd Qu.:29.50 3rd Qu.:11.5
## Max. :98.00 Max. :92.0
## NA's :48 NA's :48
## hemodialysis_years TX_order ebv cmv
## Min. : 0.000 Min. :1.000 Length:186 Length:186
## 1st Qu.: 1.400 1st Qu.:1.000 Class :character Class :character
## Median : 2.300 Median :1.000 Mode :character Mode :character
## Mean : 3.047 Mean :1.109
## 3rd Qu.: 4.000 3rd Qu.:1.000
## Max. :11.900 Max. :4.000
## NA's :49 NA's :48
## CIT_hours MT_minutes induction_therapy tacrolimus
## Min. : 0.0 Min. : 7.00 Length:186 Min. :0.0000
## 1st Qu.:12.0 1st Qu.:19.00 Class :character 1st Qu.:1.0000
## Median :14.0 Median :24.00 Mode :character Median :1.0000
## Mean :13.7 Mean :25.37 Mean :0.9783
## 3rd Qu.:16.0 3rd Qu.:29.75 3rd Qu.:1.0000
## Max. :24.0 Max. :57.00 Max. :1.0000
## NA's :48 NA's :48 NA's :48
## MMF mTOR egfr_mean creatinine_mean
## Min. :0.0000 Min. :0.00000 Min. :0.1367 Min. : 53.0
## 1st Qu.:1.0000 1st Qu.:0.00000 1st Qu.:0.4558 1st Qu.: 119.0
## Median :1.0000 Median :0.00000 Median :0.6637 Median : 254.9
## Mean :0.9275 Mean :0.02174 Mean :0.7720 Mean : 277.7
## 3rd Qu.:1.0000 3rd Qu.:0.00000 3rd Qu.:0.9631 3rd Qu.: 357.2
## Max. :1.0000 Max. :1.00000 Max. :1.8800 Max. :1172.8
## NA's :48 NA's :48
## il1_a_mean il1_b_mean il1_ra_mean il18_mean
## Min. : 0.1751 Min. : 0.3786 Min. : 1.355 Min. : 1.501
## 1st Qu.: 1.4411 1st Qu.: 8.0770 1st Qu.: 667.896 1st Qu.: 278.570
## Median : 3.7232 Median :12.7858 Median : 1147.610 Median : 431.502
## Mean : 5.0087 Mean :12.2006 Mean : 1624.957 Mean : 653.351
## 3rd Qu.: 6.5541 3rd Qu.:14.7499 3rd Qu.: 2031.337 3rd Qu.: 782.363
## Max. :46.6601 Max. :48.9390 Max. :12548.295 Max. :5015.730
##
## il18_bp_mean il18_free_mean il36_b_mean
## Min. : 1946 Min. : 77.51 Min. : 0.0258
## 1st Qu.: 2990 1st Qu.: 206.94 1st Qu.: 2.4952
## Median : 3438 Median : 317.76 Median : 3.7044
## Mean : 3587 Mean : 456.12 Mean : 5.0883
## 3rd Qu.: 4009 3rd Qu.: 565.57 3rd Qu.: 5.7242
## Max. :11804 Max. :3279.54 Max. :53.5479
## NA's :17 NA's :17
data_long <- read.table('data/long_data.txt')
data_long <- data_long%>% select(-Patient)
summary(data_long)
## id time_point group id_obs
## Length:596 Length:596 Length:596 Length:596
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## il1_a_value il1_b_value il1_ra_value il18_value
## Min. : 0.1751 Min. : 0.3786 Min. : 1.355 Min. : 1.501
## 1st Qu.: 2.2933 1st Qu.: 11.1629 1st Qu.: 665.718 1st Qu.: 279.288
## Median : 3.7308 Median : 12.7285 Median : 1032.425 Median : 448.818
## Mean : 5.5595 Mean : 13.6417 Mean : 1767.633 Mean : 715.676
## 3rd Qu.: 6.3253 3rd Qu.: 14.3496 3rd Qu.: 1745.067 3rd Qu.: 693.103
## Max. :62.5197 Max. :154.9024 Max. :24279.375 Max. :9421.688
## NA's :96 NA's :96 NA's :96 NA's :96
## il18_bp_value il18_bp_ratio il18_free il36_b_value
## Min. : 82.47 Min. : 3.765 Min. : 37.62 Min. : 0.0258
## 1st Qu.: 2976.00 1st Qu.: 31.140 1st Qu.: 193.99 1st Qu.: 2.3362
## Median : 3543.00 Median : 50.540 Median : 317.73 Median : 3.6714
## Mean : 3708.59 Mean : 87.717 Mean : 486.13 Mean : 5.5442
## 3rd Qu.: 4168.00 3rd Qu.: 81.141 3rd Qu.: 476.33 3rd Qu.: 5.5809
## Max. :30834.43 Max. :2967.066 Max. :6155.32 Max. :90.0555
## NA's :19 NA's :115 NA's :115 NA's :96
## creatinine egfr
## Min. : 53.0 Min. :0.0300
## 1st Qu.: 114.9 1st Qu.:0.1600
## Median : 169.7 Median :0.5800
## Mean : 315.9 Mean :0.6238
## 3rd Qu.: 480.1 3rd Qu.:0.9300
## Max. :1778.9 Max. :1.8800
## 2 Summary statistics
Summary table for all relevant variables
2.1 Healthy included
Open code
sumtab1 <- run(
patient_table %>% select(-id, -rejection_afterTX) %>%
tbl_summary( by='group',
type = list(hemodialysis_years ~ 'continuous',
N_mismatch ~ 'continuous')
) %>%
modify_caption('Subjects characteristics according to group (healthy vs. with acute rejection vs. no acute rejection') %>%
add_p(),
path = 'gitignore/run/sumtab1', reuse=TRUE)
sumtab1| Characteristic | acute_rejection, N = 441 | healthy, N = 481 | no_rejection, N = 941 | p-value2 |
|---|---|---|---|---|
| receiver_sex | 0.075 | |||
| female | 16 (36%) | 26 (54%) | 33 (35%) | |
| male | 28 (64%) | 22 (46%) | 61 (65%) | |
| receiver_age | 54 (44, 62) | NA (NA, NA) | 58 (46, 66) | 0.14 |
| Unknown | 0 | 48 | 0 | |
| donor_sex | >0.9 | |||
| female | 18 (41%) | 0 (NA%) | 37 (39%) | |
| male | 26 (59%) | 0 (NA%) | 57 (61%) | |
| Unknown | 0 | 48 | 0 | |
| donor_age | 51 (42, 58) | NA (NA, NA) | 57 (46, 66) | 0.066 |
| Unknown | 0 | 48 | 0 | |
| N_mismatch | 3 (2, 4) | NA (NA, NA) | 3 (3, 4) | >0.9 |
| Unknown | 0 | 48 | 0 | |
| dsa | 0.002 | |||
| neg | 27 (61%) | 0 (NA%) | 81 (86%) | |
| pos | 17 (39%) | 0 (NA%) | 13 (14%) | |
| Unknown | 0 | 48 | 0 | |
| anti_hla | 0.032 | |||
| N/A | 0 (0%) | 0 (NA%) | 1 (1.1%) | |
| neg | 27 (61%) | 0 (NA%) | 75 (80%) | |
| pos | 17 (39%) | 0 (NA%) | 18 (19%) | |
| Unknown | 0 | 48 | 0 | |
| PRA_max | 10 (3, 29) | NA (NA, NA) | 10 (2, 30) | 0.9 |
| Unknown | 0 | 48 | 0 | |
| PRA_actual | 2 (0, 6) | NA (NA, NA) | 2 (0, 20) | 0.3 |
| Unknown | 0 | 48 | 0 | |
| hemodialysis_years | 2.30 (1.40, 3.90) | NA (NA, NA) | 2.30 (1.30, 4.00) | >0.9 |
| Unknown | 0 | 48 | 1 | |
| TX_order | <0.001 | |||
| 1 | 32 (73%) | 0 (NA%) | 93 (99%) | |
| 2 | 11 (25%) | 0 (NA%) | 1 (1.1%) | |
| 4 | 1 (2.3%) | 0 (NA%) | 0 (0%) | |
| Unknown | 0 | 48 | 0 | |
| ebv | >0.9 | |||
| neg | 0 (0%) | 0 (NA%) | 2 (2.8%) | |
| pos | 35 (100%) | 0 (NA%) | 69 (97%) | |
| Unknown | 9 | 48 | 23 | |
| cmv | 0.2 | |||
| N/A | 0 (0%) | 0 (NA%) | 1 (1.1%) | |
| neg | 10 (23%) | 0 (NA%) | 11 (12%) | |
| pos | 34 (77%) | 0 (NA%) | 82 (87%) | |
| Unknown | 0 | 48 | 0 | |
| CIT_hours | 15 (12, 18) | NA (NA, NA) | 14 (12, 16) | 0.093 |
| Unknown | 0 | 48 | 0 | |
| MT_minutes | 22 (18, 27) | NA (NA, NA) | 25 (20, 31) | 0.051 |
| Unknown | 0 | 48 | 0 | |
| induction_therapy | 0.004 | |||
| ATG | 14 (32%) | 0 (NA%) | 57 (61%) | |
| ATG,PE,IVIG | 10 (23%) | 0 (NA%) | 8 (8.5%) | |
| ATG+ Infliximab | 1 (2.3%) | 0 (NA%) | 0 (0%) | |
| ATG+Infliximab | 3 (6.8%) | 0 (NA%) | 2 (2.1%) | |
| Simulect | 16 (36%) | 0 (NA%) | 27 (29%) | |
| Unknown | 0 | 48 | 0 | |
| tacrolimus | 44 (100%) | 0 (NA%) | 91 (97%) | 0.6 |
| Unknown | 0 | 48 | 0 | |
| MMF | 39 (89%) | 0 (NA%) | 89 (95%) | 0.3 |
| Unknown | 0 | 48 | 0 | |
| mTOR | 0 (0%) | 0 (NA%) | 3 (3.2%) | 0.6 |
| Unknown | 0 | 48 | 0 | |
| egfr_mean | 0.49 (0.35, 0.65) | 1.48 (1.15, 1.63) | 0.55 (0.41, 0.72) | <0.001 |
| creatinine_mean | 318 (246, 385) | 85 (71, 97) | 307 (235, 412) | <0.001 |
| il1_a_mean | 3.6 (2.9, 4.7) | 0.2 (0.2, 0.9) | 5.4 (3.5, 10.2) | <0.001 |
| il1_b_mean | 13.4 (12.5, 14.5) | 3.8 (2.9, 4.4) | 13.7 (12.2, 16.3) | <0.001 |
| il1_ra_mean | 1,291 (880, 2,044) | 440 (321, 577) | 1,418 (1,060, 2,738) | <0.001 |
| il18_mean | 626 (389, 907) | 196 (146, 280) | 554 (371, 1,010) | <0.001 |
| il18_bp_mean | 3,389 (2,999, 3,997) | 2,813 (2,465, 3,074) | 3,649 (3,346, 4,194) | <0.001 |
| Unknown | 0 | 17 | 0 | |
| il18_free_mean | 418 (268, 608) | 145 (121, 194) | 359 (236, 647) | <0.001 |
| Unknown | 0 | 17 | 0 | |
| il36_b_mean | 3.5 (2.5, 4.2) | 1.8 (1.3, 2.5) | 5.1 (3.8, 7.2) | <0.001 |
| 1 n (%); Median (IQR) | ||||
| 2 Pearson’s Chi-squared test; Kruskal-Wallis rank sum test; Fisher’s exact test | ||||
2.2 TX patients only
Open code
sumtab2 <- run(
patient_table %>% select(-id, -rejection_afterTX) %>%
filter(group != 'healthy') %>%
mutate(group = factor(group)) %>%
tbl_summary( by='group',
type = list(hemodialysis_years ~ 'continuous',
N_mismatch ~ 'continuous')
) %>%
modify_caption('Patient characteristics according to developed acute rejection') %>%
add_p(),
path = 'gitignore/run/sumtab2', reuse=TRUE)
sumtab2| Characteristic | acute_rejection, N = 441 | no_rejection, N = 941 | p-value2 |
|---|---|---|---|
| receiver_sex | 0.9 | ||
| female | 16 (36%) | 33 (35%) | |
| male | 28 (64%) | 61 (65%) | |
| receiver_age | 54 (44, 62) | 58 (46, 66) | 0.14 |
| donor_sex | 0.9 | ||
| female | 18 (41%) | 37 (39%) | |
| male | 26 (59%) | 57 (61%) | |
| donor_age | 51 (42, 58) | 57 (46, 66) | 0.066 |
| N_mismatch | 3 (2, 4) | 3 (3, 4) | >0.9 |
| dsa | <0.001 | ||
| neg | 27 (61%) | 81 (86%) | |
| pos | 17 (39%) | 13 (14%) | |
| anti_hla | 0.032 | ||
| N/A | 0 (0%) | 1 (1.1%) | |
| neg | 27 (61%) | 75 (80%) | |
| pos | 17 (39%) | 18 (19%) | |
| PRA_max | 10 (3, 29) | 10 (2, 30) | 0.9 |
| PRA_actual | 2 (0, 6) | 2 (0, 20) | 0.3 |
| hemodialysis_years | 2.30 (1.40, 3.90) | 2.30 (1.30, 4.00) | >0.9 |
| Unknown | 0 | 1 | |
| TX_order | <0.001 | ||
| 1 | 32 (73%) | 93 (99%) | |
| 2 | 11 (25%) | 1 (1.1%) | |
| 4 | 1 (2.3%) | 0 (0%) | |
| ebv | >0.9 | ||
| neg | 0 (0%) | 2 (2.8%) | |
| pos | 35 (100%) | 69 (97%) | |
| Unknown | 9 | 23 | |
| cmv | 0.2 | ||
| N/A | 0 (0%) | 1 (1.1%) | |
| neg | 10 (23%) | 11 (12%) | |
| pos | 34 (77%) | 82 (87%) | |
| CIT_hours | 15 (12, 18) | 14 (12, 16) | 0.094 |
| MT_minutes | 22 (18, 27) | 25 (20, 31) | 0.051 |
| induction_therapy | 0.004 | ||
| ATG | 14 (32%) | 57 (61%) | |
| ATG,PE,IVIG | 10 (23%) | 8 (8.5%) | |
| ATG+ Infliximab | 1 (2.3%) | 0 (0%) | |
| ATG+Infliximab | 3 (6.8%) | 2 (2.1%) | |
| Simulect | 16 (36%) | 27 (29%) | |
| tacrolimus | 44 (100%) | 91 (97%) | 0.6 |
| MMF | 39 (89%) | 89 (95%) | 0.3 |
| mTOR | 0 (0%) | 3 (3.2%) | 0.6 |
| egfr_mean | 0.49 (0.35, 0.65) | 0.55 (0.41, 0.72) | 0.2 |
| creatinine_mean | 318 (246, 385) | 307 (235, 412) | 0.9 |
| il1_a_mean | 3.6 (2.9, 4.7) | 5.4 (3.5, 10.2) | 0.001 |
| il1_b_mean | 13.4 (12.5, 14.5) | 13.7 (12.2, 16.3) | 0.13 |
| il1_ra_mean | 1,291 (880, 2,044) | 1,418 (1,060, 2,738) | 0.10 |
| il18_mean | 626 (389, 907) | 554 (371, 1,010) | >0.9 |
| il18_bp_mean | 3,389 (2,999, 3,997) | 3,649 (3,346, 4,194) | 0.2 |
| il18_free_mean | 418 (268, 608) | 359 (236, 647) | 0.7 |
| il36_b_mean | 3.5 (2.5, 4.2) | 5.1 (3.8, 7.2) | <0.001 |
| 1 n (%); Median (IQR) | |||
| 2 Pearson’s Chi-squared test; Wilcoxon rank sum test; Fisher’s exact test | |||
2.3 Summary table for continuous variables
Open code
sumtab <- patient_table %>%
select(
receiver_age, donor_age,
N_mismatch, PRA_max, PRA_actual, CIT_hours, MT_minutes,
hemodialysis_years, group, il1_a_mean, il1_b_mean, il1_ra_mean,
il18_mean, il18_bp_mean, il18_free_mean, il36_b_mean) %>%
group_by(group) %>%
summarize(across(where(is.numeric),
list(median = ~median(., na.rm = TRUE),
mean = ~mean(., na.rm = TRUE),
IQR = ~IQR(., na.rm = TRUE),
Q25 = ~quantile(., 0.25, na.rm = TRUE),
Q75 = ~quantile(., 0.75, na.rm = TRUE),
sd = ~sd(., na.rm = TRUE)),
.names = "{.col}_{.fn}"))
kableExtra::kable(t(sumtab)) | group | acute_rejection | healthy | no_rejection |
| receiver_age_median | 54 | NA | 58 |
| receiver_age_mean | 52.15909 | NA | 55.85106 |
| receiver_age_IQR | 18.25 | NA | 19.75 |
| receiver_age_Q25 | 44.00 | NA | 46.25 |
| receiver_age_Q75 | 62.25 | NA | 66.00 |
| receiver_age_sd | 13.95498 | NA | 13.41076 |
| donor_age_median | 51.0 | NA | 56.5 |
| donor_age_mean | 49.90909 | NA | 54.15957 |
| donor_age_IQR | 16.5 | NA | 20.0 |
| donor_age_Q25 | 41.75 | NA | 46.00 |
| donor_age_Q75 | 58.25 | NA | 66.00 |
| donor_age_sd | 14.00468 | NA | 15.41630 |
| N_mismatch_median | 3 | NA | 3 |
| N_mismatch_mean | 3.250000 | NA | 3.329787 |
| N_mismatch_IQR | 2 | NA | 1 |
| N_mismatch_Q25 | 2 | NA | 3 |
| N_mismatch_Q75 | 4 | NA | 4 |
| N_mismatch_sd | 1.331637 | NA | 1.212731 |
| PRA_max_median | 10 | NA | 10 |
| PRA_max_mean | 19.40909 | NA | 18.93617 |
| PRA_max_IQR | 26.25 | NA | 27.25 |
| PRA_max_Q25 | 2.75 | NA | 2.25 |
| PRA_max_Q75 | 29.0 | NA | 29.5 |
| PRA_max_sd | 23.26322 | NA | 19.97407 |
| PRA_actual_median | 2 | NA | 2 |
| PRA_actual_mean | 10.00000 | NA | 12.35106 |
| PRA_actual_IQR | 6.0 | NA | 19.5 |
| PRA_actual_Q25 | 0 | NA | 0 |
| PRA_actual_Q75 | 6.0 | NA | 19.5 |
| PRA_actual_sd | 21.58703 | NA | 18.72191 |
| CIT_hours_median | 15 | NA | 14 |
| CIT_hours_mean | 14.31818 | NA | 13.40426 |
| CIT_hours_IQR | 6 | NA | 4 |
| CIT_hours_Q25 | 12 | NA | 12 |
| CIT_hours_Q75 | 18 | NA | 16 |
| CIT_hours_sd | 5.065950 | NA | 4.489699 |
| MT_minutes_median | 21.5 | NA | 25.0 |
| MT_minutes_mean | 23.90909 | NA | 26.05319 |
| MT_minutes_IQR | 8.50 | NA | 10.75 |
| MT_minutes_Q25 | 18.00 | NA | 20.25 |
| MT_minutes_Q75 | 26.5 | NA | 31.0 |
| MT_minutes_sd | 8.454580 | NA | 8.773576 |
| hemodialysis_years_median | 2.3 | NA | 2.3 |
| hemodialysis_years_mean | 2.972727 | NA | 3.082796 |
| hemodialysis_years_IQR | 2.5 | NA | 2.7 |
| hemodialysis_years_Q25 | 1.4 | NA | 1.3 |
| hemodialysis_years_Q75 | 3.9 | NA | 4.0 |
| hemodialysis_years_sd | 2.271396 | NA | 2.491273 |
| il1_a_mean_median | 3.6138431 | 0.1751254 | 5.4099081 |
| il1_a_mean_mean | 4.519618 | 1.790396 | 6.880929 |
| il1_a_mean_IQR | 1.8399964 | 0.7684427 | 6.6405988 |
| il1_a_mean_Q25 | 2.8678507 | 0.1751254 | 3.5186279 |
| il1_a_mean_Q75 | 4.7078471 | 0.9435681 | 10.1592267 |
| il1_a_mean_sd | 3.702535 | 6.758296 | 5.079430 |
| il1_b_mean_median | 13.379882 | 3.800578 | 13.726654 |
| il1_b_mean_mean | 13.614621 | 4.568534 | 15.435899 |
| il1_b_mean_IQR | 2.016068 | 1.490545 | 4.087478 |
| il1_b_mean_Q25 | 12.456512 | 2.865737 | 12.211393 |
| il1_b_mean_Q75 | 14.472580 | 4.356282 | 16.298871 |
| il1_b_mean_sd | 1.876318 | 6.136904 | 5.552799 |
| il1_ra_mean_median | 1291.0082 | 439.8988 | 1417.6213 |
| il1_ra_mean_mean | 1607.137 | 531.280 | 2191.771 |
| il1_ra_mean_IQR | 1163.9840 | 255.8475 | 1677.9097 |
| il1_ra_mean_Q25 | 879.6144 | 320.7617 | 1059.8544 |
| il1_ra_mean_Q75 | 2043.5983 | 576.6092 | 2737.7641 |
| il1_ra_mean_sd | 1054.4684 | 464.7765 | 1883.0248 |
| il18_mean_median | 626.2038 | 196.2464 | 553.9956 |
| il18_mean_mean | 706.7674 | 220.7274 | 849.2618 |
| il18_mean_IQR | 518.0835 | 133.1681 | 639.2590 |
| il18_mean_Q25 | 388.6040 | 146.4606 | 370.9122 |
| il18_mean_Q75 | 906.6875 | 279.6288 | 1010.1712 |
| il18_mean_sd | 413.15414 | 99.47696 | 778.91956 |
| il18_bp_mean_median | 3389.167 | 2813.000 | 3649.250 |
| il18_bp_mean_mean | 3666.610 | 2854.484 | 3791.695 |
| il18_bp_mean_IQR | 998.5208 | 609.0000 | 847.5291 |
| il18_bp_mean_Q25 | 2998.875 | 2465.000 | 3346.188 |
| il18_bp_mean_Q75 | 3997.396 | 3074.000 | 4193.717 |
| il18_bp_mean_sd | 889.4910 | 514.8316 | 1053.8108 |
| il18_free_mean_median | 417.7738 | 144.8761 | 358.7065 |
| il18_free_mean_mean | 478.2362 | 164.6453 | 541.8925 |
| il18_free_mean_IQR | 339.20872 | 72.47453 | 411.11565 |
| il18_free_mean_Q25 | 268.4762 | 121.4914 | 235.7412 |
| il18_free_mean_Q75 | 607.6849 | 193.9659 | 646.8569 |
| il18_free_mean_sd | 286.58937 | 69.54324 | 499.73260 |
| il36_b_mean_median | 3.516403 | 1.823203 | 5.079123 |
| il36_b_mean_mean | 5.616041 | 2.465567 | 6.180597 |
| il36_b_mean_IQR | 1.639928 | 1.250801 | 3.346720 |
| il36_b_mean_Q25 | 2.543047 | 1.279165 | 3.806194 |
| il36_b_mean_Q75 | 4.182975 | 2.529966 | 7.152914 |
| il36_b_mean_sd | 8.701474 | 3.398531 | 4.137029 |
3 Cytokine levels changes over time and across groups
Creating a new variable, group_sep, that combines group and time_point into a single variable. Re-defining group and time_point to be factors
Open code
data_long <- data_long %>%
mutate(group_sep = factor(interaction(group, time_point)),
group = factor(group, levels =
c("healthy",
"no_rejection",
"acute_rejection")),
time_point = factor(time_point))3.1 Il1_a
3.1.1 Mixed model
Open code
## filtering data
model_data <- data_long %>%
filter(!is.na(il1_a_value))
## model fit
model_il1_a <- lmer(log10(il1_a_value) ~
group_sep +
(1|id), data = model_data)
## diagnostic plot of model
plot(model_il1_a)
Open code
## anova result od model
Anova(model_il1_a, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il1_a_value)
## F Df Df.res Pr(>F)
## group_sep 15.739 9 417.45 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## model summary
summary(model_il1_a)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il1_a_value) ~ group_sep + (1 | id)
## Data: model_data
##
## REML criterion at convergence: 807.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5688 -0.3288 0.0750 0.5766 3.4928
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.05345 0.2312
## Residual 0.23773 0.4876
## Number of obs: 500, groups: id, 186
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.54459 0.08602 6.331
## group_sephealthy.000 -0.96152 0.11604 -8.286
## group_sepno_rejection.000 0.01385 0.10245 0.135
## group_sepacute_rejection.007 -0.08337 0.10815 -0.771
## group_sepno_rejection.007 -0.09297 0.10245 -0.907
## group_sepacute_rejection.090 -0.02344 0.10975 -0.214
## group_sepno_rejection.090 0.19470 0.12090 1.610
## group_sepacute_rejection.365 -0.11856 0.11488 -1.032
## group_sepno_rejection.365 0.13162 0.12090 1.089
## group_sepacute_rejection.biopsy -0.06165 0.12267 -0.503
##
## Correlation of Fixed Effects:
## (Intr) g_.000 g__.000 grp_spc_.007 grp_spn_.007 grp_spc_.090
## grp_sph.000 -0.741
## grp_sp_.000 -0.840 0.622
## grp_spc_.007 -0.663 0.491 0.557
## grp_spn_.007 -0.840 0.622 0.759 0.557
## grp_spc_.090 -0.656 0.486 0.551 0.520 0.551
## grp_spn_.090 -0.711 0.527 0.643 0.472 0.643 0.467
## grp_spc_.365 -0.628 0.466 0.528 0.498 0.528 0.492
## grp_spn_.365 -0.711 0.527 0.643 0.472 0.643 0.467
## grp_spct_r. -0.586 0.434 0.492 0.464 0.492 0.459
## grp_spn_.090 grp_spc_.365 grp_spn_.365
## grp_sph.000
## grp_sp_.000
## grp_spc_.007
## grp_spn_.007
## grp_spc_.090
## grp_spn_.090
## grp_spc_.365 0.447
## grp_spn_.365 0.583 0.447
## grp_spct_r. 0.417 0.435 0.4173.1.2 Pairwise comparison
3.1.2.1 List of contrasts
List of contrast names and left and right coordinates on the X axis
Open code
contrast <- c(
## healthy vs. others
"healthy vs. no rejection D_0",
"healthy vs. no rejection D_7",
"healthy vs. no rejection D_90",
"healthy vs. no rejection D_365",
"healthy vs. acute rejection D_0",
"healthy vs. acute rejection D_7",
"healthy vs. acute rejection D_90",
"healthy vs. acute rejection D_365",
"healthy vs. acute rejection D_biopsy",
## between-subjects comparisons
"no rejection D_0 vs. acute rejection D_0",
"no rejection D_7 vs. acute rejection D_7",
"no rejection D_90 vs. acute rejection D_90",
"no rejection D_365 vs. acute rejection D_365",
"no rejection D_0 vs. acute rejection D_biopsy",
"no rejection D_7 vs. acute rejection D_biopsy",
"no rejection D_90 vs. acute rejection D_biopsy",
"no rejection D_365 vs. acute rejection D_biopsy",
## within-subjects comparisons
### no rejection
"no rejection D_0 vs. no rejection D_7",
"no rejection D_0 vs. no rejection D_90",
"no rejection D_0 vs. no rejection D_365",
"no rejection D_7 vs. no rejection D_90",
"no rejection D_7 vs. no rejection D_365",
"no rejection D_90 vs. no rejection D_365",
### acute rejection
"acute rejection D_0 vs. acute rejection D_7",
"acute rejection D_0 vs. acute rejection D_90",
"acute rejection D_0 vs. acute rejection D_365",
"acute rejection D_7 vs. acute rejection D_90",
"acute rejection D_7 vs. acute rejection D_365",
"acute rejection D_90 vs. acute rejection D_365",
"acute rejection D_0 vs. acute rejection D_biopsy",
"acute rejection D_7 vs. acute rejection D_biopsy",
"acute rejection D_90 vs. acute rejection D_biopsy",
"acute rejection D_365 vs. acute rejection D_biopsy"
)
xcoord_1 <- c(rep(6, 9), seq(0.8, 3.8, 1), seq(0.8, 3.8, 1),
rep(0.8, 3), rep(1.8, 2), 2.8,
rep(1.2, 3), rep(2.2, 2), 3.2,
seq(1.2, 4.2, 1))
xcoord_2 <- c(seq(0.8, 3.8, 1), seq(1.2, 4.2, 1),
5, seq(1.2, 4.2, 1), rep(5, 4),
seq(1.8, 3.8, 1), 2.8, 3.8, 3.8,
2.2, 3.2, 4.2, 3.2, 4.2, 4.2,
rep(5, 4))
anot_coord <- (xcoord_1+xcoord_2)/2
pair_table <- data.frame(contrast, xcoord_1, xcoord_2,
anot_coord)
pair_table <- pair_table %>%
mutate(term1 = if_else(xcoord_1 == 6, 'healthy.000', 'other')) %>%
mutate(term1 = if_else(xcoord_1 == 0.8, 'no_rejection.000', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 1.8, 'no_rejection.007', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 2.8, 'no_rejection.090', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 3.8, 'no_rejection.365', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 5, 'acute_rejection.biopsy', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 1.2, 'acute_rejection.000', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 2.2, 'acute_rejection.007', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 3.2, 'acute_rejection.090', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 4.2, 'acute_rejection.365', term1)) %>%
mutate(term2 = if_else(xcoord_2 == 6, 'healthy.000', 'other')) %>%
mutate(term2 = if_else(xcoord_2 == 0.8, 'no_rejection.000', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 1.8, 'no_rejection.007', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 2.8, 'no_rejection.090', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 3.8, 'no_rejection.365', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 5, 'acute_rejection.biopsy', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 1.2, 'acute_rejection.000', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 2.2, 'acute_rejection.007', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 3.2, 'acute_rejection.090', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 4.2, 'acute_rejection.365', term2))
pair_table$p <- NA3.1.2.2 Between-subjects comparisons
Open code
set.seed(16)
for (i in 1:17){
dat_pair <- model_data %>%
filter(group_sep == paste(pair_table$term1[i])|
group_sep == paste(pair_table$term2[i]))
g1 <- dat_pair %>% filter(group_sep == paste(pair_table$term1[i]))
g2 <- dat_pair %>% filter(group_sep == paste(pair_table$term2[i]))
pair_table$p[i] <- wilcox.test(g1$il1_a_value, g2$il1_a_value, paired=FALSE)$p.value
}3.1.2.3 Within-subject comparisons
Open code
set.seed(16)
for (i in 18:33){
dat_pair <- model_data %>%
filter(group_sep == paste(pair_table$term1[i]) |
group_sep == paste(pair_table$term2[i])) %>%
select(id, group_sep, il1_a_value) %>%
pivot_wider(names_from = group_sep, values_from = il1_a_value) %>%
filter(complete.cases(.)) %>% data.frame()
g1 <- dat_pair[,2]
g2 <- dat_pair[,3]
pair_table$p[i] <- wilcox.test(c(g1), c(g2), paired=TRUE)$p.value
}
pair_table <- pair_table[-c(14:15, 19:20, 22, 25:26, 28, 30:33),] %>% data.frame()
pair_table$cole <- c(rep("black", 15), rep("blue3", 3), rep("red4", 3))
pair_table$fdr <- p.adjust(pair_table$p, method = 'BH')3.1.2.4 Resulting table
Open code
## the whole `pair table`
pair_table
## contrast xcoord_1 xcoord_2 anot_coord
## 1 healthy vs. no rejection D_0 6.0 0.8 3.4
## 2 healthy vs. no rejection D_7 6.0 1.8 3.9
## 3 healthy vs. no rejection D_90 6.0 2.8 4.4
## 4 healthy vs. no rejection D_365 6.0 3.8 4.9
## 5 healthy vs. acute rejection D_0 6.0 1.2 3.6
## 6 healthy vs. acute rejection D_7 6.0 2.2 4.1
## 7 healthy vs. acute rejection D_90 6.0 3.2 4.6
## 8 healthy vs. acute rejection D_365 6.0 4.2 5.1
## 9 healthy vs. acute rejection D_biopsy 6.0 5.0 5.5
## 10 no rejection D_0 vs. acute rejection D_0 0.8 1.2 1.0
## 11 no rejection D_7 vs. acute rejection D_7 1.8 2.2 2.0
## 12 no rejection D_90 vs. acute rejection D_90 2.8 3.2 3.0
## 13 no rejection D_365 vs. acute rejection D_365 3.8 4.2 4.0
## 16 no rejection D_90 vs. acute rejection D_biopsy 2.8 5.0 3.9
## 17 no rejection D_365 vs. acute rejection D_biopsy 3.8 5.0 4.4
## 18 no rejection D_0 vs. no rejection D_7 0.8 1.8 1.3
## 21 no rejection D_7 vs. no rejection D_90 1.8 2.8 2.3
## 23 no rejection D_90 vs. no rejection D_365 2.8 3.8 3.3
## 24 acute rejection D_0 vs. acute rejection D_7 1.2 2.2 1.7
## 27 acute rejection D_7 vs. acute rejection D_90 2.2 3.2 2.7
## 29 acute rejection D_90 vs. acute rejection D_365 3.2 4.2 3.7
## term1 term2 p cole fdr
## 1 healthy.000 no_rejection.000 6.973501e-11 black 2.928870e-10
## 2 healthy.000 no_rejection.007 5.227359e-11 black 2.744363e-10
## 3 healthy.000 no_rejection.090 1.475037e-13 black 3.097579e-12
## 4 healthy.000 no_rejection.365 7.054097e-13 black 7.406802e-12
## 5 healthy.000 acute_rejection.000 5.467856e-10 black 1.913750e-09
## 6 healthy.000 acute_rejection.007 2.952645e-09 black 7.750693e-09
## 7 healthy.000 acute_rejection.090 3.797214e-11 black 2.658050e-10
## 8 healthy.000 acute_rejection.365 2.847456e-09 black 7.750693e-09
## 9 healthy.000 acute_rejection.biopsy 4.447710e-09 black 1.037799e-08
## 10 no_rejection.000 acute_rejection.000 3.608654e-01 black 4.300071e-01
## 11 no_rejection.007 acute_rejection.007 3.455284e-01 black 4.300071e-01
## 12 no_rejection.090 acute_rejection.090 6.924033e-06 black 1.321861e-05
## 13 no_rejection.365 acute_rejection.365 2.968675e-06 black 6.234218e-06
## 16 no_rejection.090 acute_rejection.biopsy 1.026038e-04 black 1.795567e-04
## 17 no_rejection.365 acute_rejection.biopsy 7.185841e-04 black 1.160790e-03
## 18 no_rejection.000 no_rejection.007 7.523005e-03 blue3 1.128451e-02
## 21 no_rejection.007 no_rejection.090 6.940287e-01 blue3 7.670843e-01
## 23 no_rejection.090 no_rejection.365 3.685775e-01 blue3 4.300071e-01
## 24 acute_rejection.000 acute_rejection.007 8.076073e-01 red4 8.076073e-01
## 27 acute_rejection.007 acute_rejection.090 7.446044e-01 red4 7.818346e-01
## 29 acute_rejection.090 acute_rejection.365 1.794364e-01 red4 2.512110e-01
## only significant differences
pair_table %>%
filter(p<0.05) %>%
select(contrast, p, fdr) %>%
mutate(star = if_else(p<0.05, '*', '')) %>%
mutate(star = if_else(p<0.01, '**', star)) %>%
mutate(star = if_else(p<0.001, '***', star)) %>%
data.frame()
## contrast p fdr
## 1 healthy vs. no rejection D_0 6.973501e-11 2.928870e-10
## 2 healthy vs. no rejection D_7 5.227359e-11 2.744363e-10
## 3 healthy vs. no rejection D_90 1.475037e-13 3.097579e-12
## 4 healthy vs. no rejection D_365 7.054097e-13 7.406802e-12
## 5 healthy vs. acute rejection D_0 5.467856e-10 1.913750e-09
## 6 healthy vs. acute rejection D_7 2.952645e-09 7.750693e-09
## 7 healthy vs. acute rejection D_90 3.797214e-11 2.658050e-10
## 8 healthy vs. acute rejection D_365 2.847456e-09 7.750693e-09
## 9 healthy vs. acute rejection D_biopsy 4.447710e-09 1.037799e-08
## 10 no rejection D_90 vs. acute rejection D_90 6.924033e-06 1.321861e-05
## 11 no rejection D_365 vs. acute rejection D_365 2.968675e-06 6.234218e-06
## 12 no rejection D_90 vs. acute rejection D_biopsy 1.026038e-04 1.795567e-04
## 13 no rejection D_365 vs. acute rejection D_biopsy 7.185841e-04 1.160790e-03
## 14 no rejection D_0 vs. no rejection D_7 7.523005e-03 1.128451e-02
## star
## 1 ***
## 2 ***
## 3 ***
## 4 ***
## 5 ***
## 6 ***
## 7 ***
## 8 ***
## 9 ***
## 10 ***
## 11 ***
## 12 ***
## 13 ***
## 14 **
pair_rest <- pair_table %>%
filter(p<0.05,
term1 != 'healthy.000') %>%
mutate(star = if_else(p<0.05, '*', '')) %>%
mutate(star = if_else(p<0.01, '**', star)) %>%
mutate(star = if_else(p<0.001, '***', star)) %>%
data.frame()3.1.3 Visualization
Open code
## plot
plotos <- data_long %>%
filter(!is.na(il1_a_value)) %>%
mutate(group = factor(group, levels =
c("no_rejection",
"acute_rejection",
"healthy")),
time_point = factor(time_point)) %>%
mutate(time_point = if_else(group == 'healthy', 'healthy', time_point)) %>%
mutate(group_sep = interaction(group, time_point)) %>%
## plotting
ggplot(aes(x=time_point, y=il1_a_value)) +
geom_boxplot(aes(fill=group, group = group_sep), outlier.shape = NA,
size=0.5,position = position_dodge(width = 0.84), color='black') +
labs(x = "Time post TX", y = expression("IL-1" * alpha * " (pg/ml)")) +
scale_x_discrete(labels = c("Before TX", "7 days", "3 months", "365 days",
"Biopsy", "Healthy")) +
theme(axis.text = element_text(size = 11),
axis.title = element_text(size = 12),
legend.text = element_text(size = 11),
legend.title = element_text(size = 12)) +
scale_fill_manual(values = c("no_rejection" = "skyblue2",
"acute_rejection" = "coral",
"healthy" = "darkolivegreen3"),
labels = c("no rejection", "acute rejection", "healthy")) +
coord_cartesian(ylim = c(0, 30)) +
annotate("text", x = 6, y = 5, label = "#", size = 7)
pair_rest$y <- c(12, 11, 16, 13,30)
plotos <- plotos + geom_segment(data = pair_rest,
aes(x = xcoord_1, xend = xcoord_2, y = y, yend = y),
color = pair_rest$cole,
linetype = "solid", size = 0.4)
plot_il1_a <- plotos + annotate("text", x = pair_rest$anot_coord, y = pair_rest$y+0.5,
label = pair_rest$star, size = 7, color=pair_rest$cole)
plot_il1_a
3.2 Il1_b
3.2.1 Mixed model
Open code
## filtering data
model_data <- data_long %>%
filter(!is.na(il1_b_value))
## model fit
model_il1_b <- lmer(log10(il1_b_value) ~
group_sep +
(1|id), data = model_data)
## diagnostic plot of model
plot(model_il1_b)
Open code
## anova result od model
Anova(model_il1_b, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il1_b_value)
## F Df Df.res Pr(>F)
## group_sep 83.17 9 417.24 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## model summary
summary(model_il1_b)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il1_b_value) ~ group_sep + (1 | id)
## Data: model_data
##
## REML criterion at convergence: -509.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.2908 -0.5140 -0.0895 0.2629 7.1938
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.003725 0.06103
## Residual 0.016117 0.12695
## Number of obs: 500, groups: id, 186
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.16016 0.02245 51.675
## group_sephealthy.000 -0.59874 0.03029 -19.768
## group_sepno_rejection.000 0.04741 0.02674 1.773
## group_sepacute_rejection.007 -0.05044 0.02816 -1.791
## group_sepno_rejection.007 -0.01361 0.02674 -0.509
## group_sepacute_rejection.090 -0.04086 0.02858 -1.430
## group_sepno_rejection.090 -0.07554 0.03155 -2.395
## group_sepacute_rejection.365 -0.05930 0.02992 -1.982
## group_sepno_rejection.365 -0.09149 0.03155 -2.900
## group_sepacute_rejection.biopsy -0.04640 0.03195 -1.453
##
## Correlation of Fixed Effects:
## (Intr) g_.000 g__.000 grp_spc_.007 grp_spn_.007 grp_spc_.090
## grp_sph.000 -0.741
## grp_sp_.000 -0.840 0.622
## grp_spc_.007 -0.661 0.490 0.555
## grp_spn_.007 -0.840 0.622 0.760 0.555
## grp_spc_.090 -0.654 0.485 0.549 0.520 0.549
## grp_spn_.090 -0.712 0.528 0.644 0.471 0.644 0.466
## grp_spc_.365 -0.627 0.465 0.526 0.498 0.526 0.492
## grp_spn_.365 -0.712 0.528 0.644 0.471 0.644 0.466
## grp_spct_r. -0.584 0.433 0.490 0.464 0.490 0.459
## grp_spn_.090 grp_spc_.365 grp_spn_.365
## grp_sph.000
## grp_sp_.000
## grp_spc_.007
## grp_spn_.007
## grp_spc_.090
## grp_spn_.090
## grp_spc_.365 0.446
## grp_spn_.365 0.585 0.446
## grp_spct_r. 0.416 0.435 0.4163.2.2 Pairwise comparison
3.2.2.1 List of contrasts
List of contrast names and left and right coordinates on the X axis
Open code
contrast <- c(
## healthy vs. others
"healthy vs. no rejection D_0",
"healthy vs. no rejection D_7",
"healthy vs. no rejection D_90",
"healthy vs. no rejection D_365",
"healthy vs. acute rejection D_0",
"healthy vs. acute rejection D_7",
"healthy vs. acute rejection D_90",
"healthy vs. acute rejection D_365",
"healthy vs. acute rejection D_biopsy",
## between-subjects comparisons
"no rejection D_0 vs. acute rejection D_0",
"no rejection D_7 vs. acute rejection D_7",
"no rejection D_90 vs. acute rejection D_90",
"no rejection D_365 vs. acute rejection D_365",
"no rejection D_0 vs. acute rejection D_biopsy",
"no rejection D_7 vs. acute rejection D_biopsy",
"no rejection D_90 vs. acute rejection D_biopsy",
"no rejection D_365 vs. acute rejection D_biopsy",
## within-subjects comparisons
### no rejection
"no rejection D_0 vs. no rejection D_7",
"no rejection D_0 vs. no rejection D_90",
"no rejection D_0 vs. no rejection D_365",
"no rejection D_7 vs. no rejection D_90",
"no rejection D_7 vs. no rejection D_365",
"no rejection D_90 vs. no rejection D_365",
### acute rejection
"acute rejection D_0 vs. acute rejection D_7",
"acute rejection D_0 vs. acute rejection D_90",
"acute rejection D_0 vs. acute rejection D_365",
"acute rejection D_7 vs. acute rejection D_90",
"acute rejection D_7 vs. acute rejection D_365",
"acute rejection D_90 vs. acute rejection D_365",
"acute rejection D_0 vs. acute rejection D_biopsy",
"acute rejection D_7 vs. acute rejection D_biopsy",
"acute rejection D_90 vs. acute rejection D_biopsy",
"acute rejection D_365 vs. acute rejection D_biopsy"
)
xcoord_1 <- c(rep(6, 9), seq(0.8, 3.8, 1), seq(0.8, 3.8, 1),
rep(0.8, 3), rep(1.8, 2), 2.8,
rep(1.2, 3), rep(2.2, 2), 3.2,
seq(1.2, 4.2, 1))
xcoord_2 <- c(seq(0.8, 3.8, 1), seq(1.2, 4.2, 1),
5, seq(1.2, 4.2, 1), rep(5, 4),
seq(1.8, 3.8, 1), 2.8, 3.8, 3.8,
2.2, 3.2, 4.2, 3.2, 4.2, 4.2,
rep(5, 4))
anot_coord <- (xcoord_1+xcoord_2)/2
pair_table <- data.frame(contrast, xcoord_1, xcoord_2,
anot_coord)
pair_table <- pair_table %>%
mutate(term1 = if_else(xcoord_1 == 6, 'healthy.000', 'other')) %>%
mutate(term1 = if_else(xcoord_1 == 0.8, 'no_rejection.000', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 1.8, 'no_rejection.007', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 2.8, 'no_rejection.090', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 3.8, 'no_rejection.365', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 5, 'acute_rejection.biopsy', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 1.2, 'acute_rejection.000', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 2.2, 'acute_rejection.007', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 3.2, 'acute_rejection.090', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 4.2, 'acute_rejection.365', term1)) %>%
mutate(term2 = if_else(xcoord_2 == 6, 'healthy.000', 'other')) %>%
mutate(term2 = if_else(xcoord_2 == 0.8, 'no_rejection.000', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 1.8, 'no_rejection.007', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 2.8, 'no_rejection.090', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 3.8, 'no_rejection.365', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 5, 'acute_rejection.biopsy', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 1.2, 'acute_rejection.000', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 2.2, 'acute_rejection.007', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 3.2, 'acute_rejection.090', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 4.2, 'acute_rejection.365', term2))
pair_table$p <- NA3.2.2.2 Between-subjects comparisons
Open code
set.seed(16)
for (i in 1:17){
dat_pair <- model_data %>%
filter(group_sep == paste(pair_table$term1[i])|
group_sep == paste(pair_table$term2[i]))
g1 <- dat_pair %>% filter(group_sep == paste(pair_table$term1[i]))
g2 <- dat_pair %>% filter(group_sep == paste(pair_table$term2[i]))
pair_table$p[i] <- wilcox.test(g1$il1_b_value, g2$il1_b_value, paired=FALSE)$p.value
}3.2.2.3 Within-subject comparisons
Open code
set.seed(16)
for (i in 18:33){
dat_pair <- model_data %>%
filter(group_sep == paste(pair_table$term1[i]) |
group_sep == paste(pair_table$term2[i])) %>%
select(id, group_sep, il1_b_value) %>%
pivot_wider(names_from = group_sep, values_from = il1_b_value) %>%
filter(complete.cases(.)) %>% data.frame()
g1 <- dat_pair[,2]
g2 <- dat_pair[,3]
pair_table$p[i] <- wilcox.test(c(g1), c(g2), paired=TRUE)$p.value
}
pair_table <- pair_table[-c(14:15, 19:20, 22, 25:26, 28, 30:33),]
pair_table$cole <- c(rep("black", 15), rep("blue3", 3), rep("red4", 3))
pair_table$fdr <- p.adjust(pair_table$p, method = 'BH')3.2.2.4 Resulting table
Open code
## the whole `pair table`
pair_table
## contrast xcoord_1 xcoord_2 anot_coord
## 1 healthy vs. no rejection D_0 6.0 0.8 3.4
## 2 healthy vs. no rejection D_7 6.0 1.8 3.9
## 3 healthy vs. no rejection D_90 6.0 2.8 4.4
## 4 healthy vs. no rejection D_365 6.0 3.8 4.9
## 5 healthy vs. acute rejection D_0 6.0 1.2 3.6
## 6 healthy vs. acute rejection D_7 6.0 2.2 4.1
## 7 healthy vs. acute rejection D_90 6.0 3.2 4.6
## 8 healthy vs. acute rejection D_365 6.0 4.2 5.1
## 9 healthy vs. acute rejection D_biopsy 6.0 5.0 5.5
## 10 no rejection D_0 vs. acute rejection D_0 0.8 1.2 1.0
## 11 no rejection D_7 vs. acute rejection D_7 1.8 2.2 2.0
## 12 no rejection D_90 vs. acute rejection D_90 2.8 3.2 3.0
## 13 no rejection D_365 vs. acute rejection D_365 3.8 4.2 4.0
## 16 no rejection D_90 vs. acute rejection D_biopsy 2.8 5.0 3.9
## 17 no rejection D_365 vs. acute rejection D_biopsy 3.8 5.0 4.4
## 18 no rejection D_0 vs. no rejection D_7 0.8 1.8 1.3
## 21 no rejection D_7 vs. no rejection D_90 1.8 2.8 2.3
## 23 no rejection D_90 vs. no rejection D_365 2.8 3.8 3.3
## 24 acute rejection D_0 vs. acute rejection D_7 1.2 2.2 1.7
## 27 acute rejection D_7 vs. acute rejection D_90 2.2 3.2 2.7
## 29 acute rejection D_90 vs. acute rejection D_365 3.2 4.2 3.7
## term1 term2 p cole fdr
## 1 healthy.000 no_rejection.000 1.007494e-20 black 1.100496e-19
## 2 healthy.000 no_rejection.007 1.048092e-20 black 1.100496e-19
## 3 healthy.000 no_rejection.090 1.941165e-14 black 6.224985e-14
## 4 healthy.000 no_rejection.365 2.074995e-14 black 6.224985e-14
## 5 healthy.000 acute_rejection.000 1.954387e-14 black 6.224985e-14
## 6 healthy.000 acute_rejection.007 3.839545e-15 black 2.687682e-14
## 7 healthy.000 acute_rejection.090 8.457201e-15 black 4.440030e-14
## 8 healthy.000 acute_rejection.365 1.162541e-13 black 3.051669e-13
## 9 healthy.000 acute_rejection.biopsy 4.140081e-12 black 9.660189e-12
## 10 no_rejection.000 acute_rejection.000 1.656273e-01 black 2.045985e-01
## 11 no_rejection.007 acute_rejection.007 3.029547e-02 black 4.241366e-02
## 12 no_rejection.090 acute_rejection.090 1.269805e-01 black 1.666619e-01
## 13 no_rejection.365 acute_rejection.365 1.059940e-02 black 1.712211e-02
## 16 no_rejection.090 acute_rejection.biopsy 2.524270e-01 black 2.944981e-01
## 17 no_rejection.365 acute_rejection.biopsy 2.029089e-02 black 3.043634e-02
## 18 no_rejection.000 no_rejection.007 8.577426e-05 blue3 1.801260e-04
## 21 no_rejection.007 no_rejection.090 1.573119e-03 blue3 3.003227e-03
## 23 no_rejection.090 no_rejection.365 5.391972e-01 blue3 5.391972e-01
## 24 acute_rejection.000 acute_rejection.007 5.187312e-03 red4 9.077796e-03
## 27 acute_rejection.007 acute_rejection.090 5.013747e-01 red4 5.264435e-01
## 29 acute_rejection.090 acute_rejection.365 4.890153e-01 red4 5.264435e-01
## only significant differences
pair_table %>% filter(p<0.05) %>% data.frame()
## contrast xcoord_1 xcoord_2 anot_coord
## 1 healthy vs. no rejection D_0 6.0 0.8 3.4
## 2 healthy vs. no rejection D_7 6.0 1.8 3.9
## 3 healthy vs. no rejection D_90 6.0 2.8 4.4
## 4 healthy vs. no rejection D_365 6.0 3.8 4.9
## 5 healthy vs. acute rejection D_0 6.0 1.2 3.6
## 6 healthy vs. acute rejection D_7 6.0 2.2 4.1
## 7 healthy vs. acute rejection D_90 6.0 3.2 4.6
## 8 healthy vs. acute rejection D_365 6.0 4.2 5.1
## 9 healthy vs. acute rejection D_biopsy 6.0 5.0 5.5
## 10 no rejection D_7 vs. acute rejection D_7 1.8 2.2 2.0
## 11 no rejection D_365 vs. acute rejection D_365 3.8 4.2 4.0
## 12 no rejection D_365 vs. acute rejection D_biopsy 3.8 5.0 4.4
## 13 no rejection D_0 vs. no rejection D_7 0.8 1.8 1.3
## 14 no rejection D_7 vs. no rejection D_90 1.8 2.8 2.3
## 15 acute rejection D_0 vs. acute rejection D_7 1.2 2.2 1.7
## term1 term2 p cole fdr
## 1 healthy.000 no_rejection.000 1.007494e-20 black 1.100496e-19
## 2 healthy.000 no_rejection.007 1.048092e-20 black 1.100496e-19
## 3 healthy.000 no_rejection.090 1.941165e-14 black 6.224985e-14
## 4 healthy.000 no_rejection.365 2.074995e-14 black 6.224985e-14
## 5 healthy.000 acute_rejection.000 1.954387e-14 black 6.224985e-14
## 6 healthy.000 acute_rejection.007 3.839545e-15 black 2.687682e-14
## 7 healthy.000 acute_rejection.090 8.457201e-15 black 4.440030e-14
## 8 healthy.000 acute_rejection.365 1.162541e-13 black 3.051669e-13
## 9 healthy.000 acute_rejection.biopsy 4.140081e-12 black 9.660189e-12
## 10 no_rejection.007 acute_rejection.007 3.029547e-02 black 4.241366e-02
## 11 no_rejection.365 acute_rejection.365 1.059940e-02 black 1.712211e-02
## 12 no_rejection.365 acute_rejection.biopsy 2.029089e-02 black 3.043634e-02
## 13 no_rejection.000 no_rejection.007 8.577426e-05 blue3 1.801260e-04
## 14 no_rejection.007 no_rejection.090 1.573119e-03 blue3 3.003227e-03
## 15 acute_rejection.000 acute_rejection.007 5.187312e-03 red4 9.077796e-03
pair_rest <- pair_table %>%
filter(p<0.05,
term1 != 'healthy.000') %>%
mutate(star = if_else(p<0.05, '*', '')) %>%
mutate(star = if_else(p<0.01, '**', star)) %>%
mutate(star = if_else(p<0.001, '***', star)) %>%
data.frame()3.2.3 Visualization
Open code
## plot
plotos <- data_long %>%
filter(!is.na(il1_b_value)) %>%
mutate(group = factor(group, levels =
c("no_rejection",
"acute_rejection",
"healthy")),
time_point = factor(time_point)) %>%
mutate(time_point = if_else(group == 'healthy', 'healthy', time_point)) %>%
mutate(group_sep = interaction(group, time_point)) %>%
## plotting
ggplot(aes(x=time_point, y=il1_b_value)) +
geom_boxplot(aes(fill=group, group = group_sep), outlier.shape = NA,
size=0.5,position = position_dodge(width = 0.84), color='black') +
labs(x = "Time post TX", y = expression("IL-1" * beta * " (pg/ml)")) +
scale_x_discrete(labels = c("Before TX", "7 days", "3 months", "365 days",
"Biopsy", "Healthy")) +
theme(axis.text = element_text(size = 11),
axis.title = element_text(size = 12),
legend.text = element_text(size = 11),
legend.title = element_text(size = 12)) +
scale_fill_manual(values = c("no_rejection" = "skyblue2",
"acute_rejection" = "coral",
"healthy" = "darkolivegreen3"),
labels = c("no rejection", "acute rejection", "healthy")) +
coord_cartesian(ylim = c(0, 28)) +
annotate("text", x = 6, y = 10, label = "#", size = 7)
pair_rest$y <- c(19, 17, 19, 26.5, 21, 25)
plotos <- plotos + geom_segment(data = pair_rest,
aes(x = xcoord_1, xend = xcoord_2, y = y, yend = y),
color = pair_rest$cole,
linetype = "solid", size = 0.4)
plot_il1_b <- plotos + annotate("text", x = pair_rest$anot_coord, y = pair_rest$y+0.5,
label = pair_rest$star, size = 7, color=pair_rest$cole)
plot_il1_b
3.3 il1_ra
3.3.1 Mixed model
Open code
## filtering data
model_data <- data_long %>%
filter(!is.na(il1_ra_value))
## model fit
model_il1_ra <- lmer(log10(il1_ra_value) ~
group_sep +
(1|id), data = model_data)
## diagnostic plot of model
plot(model_il1_ra)
Open code
## anova result od model
Anova(model_il1_ra, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il1_ra_value)
## F Df Df.res Pr(>F)
## group_sep 17.915 9 414.32 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## model summary
summary(model_il1_ra)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il1_ra_value) ~ group_sep + (1 | id)
## Data: model_data
##
## REML criterion at convergence: 325.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.3924 -0.5919 -0.1109 0.4126 3.4509
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.02780 0.1667
## Residual 0.08388 0.2896
## Number of obs: 500, groups: id, 186
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.17638 0.05311 59.805
## group_sephealthy.000 -0.57937 0.07175 -8.075
## group_sepno_rejection.000 0.11718 0.06332 1.851
## group_sepacute_rejection.007 -0.12171 0.06428 -1.893
## group_sepno_rejection.007 -0.06739 0.06332 -1.064
## group_sepacute_rejection.090 -0.11194 0.06528 -1.715
## group_sepno_rejection.090 -0.06636 0.07428 -0.893
## group_sepacute_rejection.365 -0.28796 0.06841 -4.209
## group_sepno_rejection.365 -0.17831 0.07428 -2.400
## group_sepacute_rejection.biopsy -0.17195 0.07310 -2.352
##
## Correlation of Fixed Effects:
## (Intr) g_.000 g__.000 grp_spc_.007 grp_spn_.007 grp_spc_.090
## grp_sph.000 -0.740
## grp_sp_.000 -0.839 0.621
## grp_spc_.007 -0.639 0.473 0.536
## grp_spn_.007 -0.839 0.621 0.777 0.536
## grp_spc_.090 -0.632 0.468 0.530 0.520 0.530
## grp_spn_.090 -0.715 0.529 0.663 0.457 0.663 0.452
## grp_spc_.365 -0.605 0.448 0.508 0.498 0.508 0.492
## grp_spn_.365 -0.715 0.529 0.663 0.457 0.663 0.452
## grp_spct_r. -0.563 0.417 0.472 0.463 0.472 0.459
## grp_spn_.090 grp_spc_.365 grp_spn_.365
## grp_sph.000
## grp_sp_.000
## grp_spc_.007
## grp_spn_.007
## grp_spc_.090
## grp_spn_.090
## grp_spc_.365 0.433
## grp_spn_.365 0.610 0.433
## grp_spct_r. 0.403 0.433 0.4033.3.2 Pairwise comparison
3.3.2.1 List of contrasts
List of contrast names and left and right coordinates on the X axis
Open code
contrast <- c(
## healthy vs. others
"healthy vs. no rejection D_0",
"healthy vs. no rejection D_7",
"healthy vs. no rejection D_90",
"healthy vs. no rejection D_365",
"healthy vs. acute rejection D_0",
"healthy vs. acute rejection D_7",
"healthy vs. acute rejection D_90",
"healthy vs. acute rejection D_365",
"healthy vs. acute rejection D_biopsy",
## between-subjects comparisons
"no rejection D_0 vs. acute rejection D_0",
"no rejection D_7 vs. acute rejection D_7",
"no rejection D_90 vs. acute rejection D_90",
"no rejection D_365 vs. acute rejection D_365",
"no rejection D_0 vs. acute rejection D_biopsy",
"no rejection D_7 vs. acute rejection D_biopsy",
"no rejection D_90 vs. acute rejection D_biopsy",
"no rejection D_365 vs. acute rejection D_biopsy",
## within-subjects comparisons
### no rejection
"no rejection D_0 vs. no rejection D_7",
"no rejection D_0 vs. no rejection D_90",
"no rejection D_0 vs. no rejection D_365",
"no rejection D_7 vs. no rejection D_90",
"no rejection D_7 vs. no rejection D_365",
"no rejection D_90 vs. no rejection D_365",
### acute rejection
"acute rejection D_0 vs. acute rejection D_7",
"acute rejection D_0 vs. acute rejection D_90",
"acute rejection D_0 vs. acute rejection D_365",
"acute rejection D_7 vs. acute rejection D_90",
"acute rejection D_7 vs. acute rejection D_365",
"acute rejection D_90 vs. acute rejection D_365",
"acute rejection D_0 vs. acute rejection D_biopsy",
"acute rejection D_7 vs. acute rejection D_biopsy",
"acute rejection D_90 vs. acute rejection D_biopsy",
"acute rejection D_365 vs. acute rejection D_biopsy"
)
xcoord_1 <- c(rep(6, 9), seq(0.8, 3.8, 1), seq(0.8, 3.8, 1),
rep(0.8, 3), rep(1.8, 2), 2.8,
rep(1.2, 3), rep(2.2, 2), 3.2,
seq(1.2, 4.2, 1))
xcoord_2 <- c(seq(0.8, 3.8, 1), seq(1.2, 4.2, 1),
5, seq(1.2, 4.2, 1), rep(5, 4),
seq(1.8, 3.8, 1), 2.8, 3.8, 3.8,
2.2, 3.2, 4.2, 3.2, 4.2, 4.2,
rep(5, 4))
anot_coord <- (xcoord_1+xcoord_2)/2
pair_table <- data.frame(contrast, xcoord_1, xcoord_2,
anot_coord)
pair_table <- pair_table %>%
mutate(term1 = if_else(xcoord_1 == 6, 'healthy.000', 'other')) %>%
mutate(term1 = if_else(xcoord_1 == 0.8, 'no_rejection.000', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 1.8, 'no_rejection.007', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 2.8, 'no_rejection.090', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 3.8, 'no_rejection.365', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 5, 'acute_rejection.biopsy', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 1.2, 'acute_rejection.000', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 2.2, 'acute_rejection.007', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 3.2, 'acute_rejection.090', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 4.2, 'acute_rejection.365', term1)) %>%
mutate(term2 = if_else(xcoord_2 == 6, 'healthy.000', 'other')) %>%
mutate(term2 = if_else(xcoord_2 == 0.8, 'no_rejection.000', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 1.8, 'no_rejection.007', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 2.8, 'no_rejection.090', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 3.8, 'no_rejection.365', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 5, 'acute_rejection.biopsy', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 1.2, 'acute_rejection.000', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 2.2, 'acute_rejection.007', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 3.2, 'acute_rejection.090', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 4.2, 'acute_rejection.365', term2))
pair_table$p <- NA3.3.2.2 Between-subjects comparisons
Open code
set.seed(16)
for (i in 1:17){
dat_pair <- model_data %>%
filter(group_sep == paste(pair_table$term1[i])|
group_sep == paste(pair_table$term2[i]))
g1 <- dat_pair %>% filter(group_sep == paste(pair_table$term1[i]))
g2 <- dat_pair %>% filter(group_sep == paste(pair_table$term2[i]))
pair_table$p[i] <- wilcox.test(g1$il1_ra_value, g2$il1_ra_value, paired=FALSE)$p.value
}3.3.2.3 Within-subject comparisons
Open code
set.seed(16)
for (i in 18:33){
dat_pair <- model_data %>%
filter(group_sep == paste(pair_table$term1[i]) |
group_sep == paste(pair_table$term2[i])) %>%
select(id, group_sep, il1_ra_value) %>%
pivot_wider(names_from = group_sep, values_from = il1_ra_value) %>%
filter(complete.cases(.)) %>% data.frame()
g1 <- dat_pair[,2]
g2 <- dat_pair[,3]
pair_table$p[i] <- wilcox.test(c(g1), c(g2), paired=TRUE)$p.value
}
pair_table <- pair_table[-c(14:15, 19:20, 22, 25:26, 28, 30:33),]
pair_table$cole <- c(rep("black", 15), rep("blue3", 3), rep("red4", 3))
pair_table$fdr <- p.adjust(pair_table$p, method = 'BH')3.3.2.4 Resulting table
Open code
## the whole `pair table`
pair_table
## contrast xcoord_1 xcoord_2 anot_coord
## 1 healthy vs. no rejection D_0 6.0 0.8 3.4
## 2 healthy vs. no rejection D_7 6.0 1.8 3.9
## 3 healthy vs. no rejection D_90 6.0 2.8 4.4
## 4 healthy vs. no rejection D_365 6.0 3.8 4.9
## 5 healthy vs. acute rejection D_0 6.0 1.2 3.6
## 6 healthy vs. acute rejection D_7 6.0 2.2 4.1
## 7 healthy vs. acute rejection D_90 6.0 3.2 4.6
## 8 healthy vs. acute rejection D_365 6.0 4.2 5.1
## 9 healthy vs. acute rejection D_biopsy 6.0 5.0 5.5
## 10 no rejection D_0 vs. acute rejection D_0 0.8 1.2 1.0
## 11 no rejection D_7 vs. acute rejection D_7 1.8 2.2 2.0
## 12 no rejection D_90 vs. acute rejection D_90 2.8 3.2 3.0
## 13 no rejection D_365 vs. acute rejection D_365 3.8 4.2 4.0
## 16 no rejection D_90 vs. acute rejection D_biopsy 2.8 5.0 3.9
## 17 no rejection D_365 vs. acute rejection D_biopsy 3.8 5.0 4.4
## 18 no rejection D_0 vs. no rejection D_7 0.8 1.8 1.3
## 21 no rejection D_7 vs. no rejection D_90 1.8 2.8 2.3
## 23 no rejection D_90 vs. no rejection D_365 2.8 3.8 3.3
## 24 acute rejection D_0 vs. acute rejection D_7 1.2 2.2 1.7
## 27 acute rejection D_7 vs. acute rejection D_90 2.2 3.2 2.7
## 29 acute rejection D_90 vs. acute rejection D_365 3.2 4.2 3.7
## term1 term2 p cole fdr
## 1 healthy.000 no_rejection.000 6.278863e-18 black 1.318561e-16
## 2 healthy.000 no_rejection.007 4.243487e-16 black 4.455661e-15
## 3 healthy.000 no_rejection.090 1.424391e-11 black 9.970737e-11
## 4 healthy.000 no_rejection.365 6.321421e-09 black 1.896426e-08
## 5 healthy.000 acute_rejection.000 1.340189e-10 black 7.035992e-10
## 6 healthy.000 acute_rejection.007 2.975331e-10 black 1.249639e-09
## 7 healthy.000 acute_rejection.090 2.122261e-09 black 7.427913e-09
## 8 healthy.000 acute_rejection.365 9.383021e-06 black 2.189372e-05
## 9 healthy.000 acute_rejection.biopsy 1.355816e-07 black 3.559017e-07
## 10 no_rejection.000 acute_rejection.000 1.894461e-01 black 2.340216e-01
## 11 no_rejection.007 acute_rejection.007 6.256756e-02 black 9.581960e-02
## 12 no_rejection.090 acute_rejection.090 2.120715e-01 black 2.474167e-01
## 13 no_rejection.365 acute_rejection.365 6.387973e-02 black 9.581960e-02
## 16 no_rejection.090 acute_rejection.biopsy 1.110122e-01 black 1.554170e-01
## 17 no_rejection.365 acute_rejection.biopsy 8.746157e-01 black 8.746157e-01
## 18 no_rejection.000 no_rejection.007 6.409870e-05 blue3 1.346073e-04
## 21 no_rejection.007 no_rejection.090 5.571619e-01 blue3 6.158106e-01
## 23 no_rejection.090 no_rejection.365 2.581426e-02 blue3 4.517495e-02
## 24 acute_rejection.000 acute_rejection.007 1.422637e-01 red4 1.867211e-01
## 27 acute_rejection.007 acute_rejection.090 7.546941e-01 red4 7.924288e-01
## 29 acute_rejection.090 acute_rejection.365 2.400159e-02 red4 4.517495e-02
## only significant differences
pair_table %>%
filter(p<0.05) %>%
select(contrast, p, fdr) %>%
mutate(star = if_else(p<0.05, '*', '')) %>%
mutate(star = if_else(p<0.01, '**', star)) %>%
mutate(star = if_else(p<0.001, '***', star)) %>%
data.frame()
## contrast p fdr
## 1 healthy vs. no rejection D_0 6.278863e-18 1.318561e-16
## 2 healthy vs. no rejection D_7 4.243487e-16 4.455661e-15
## 3 healthy vs. no rejection D_90 1.424391e-11 9.970737e-11
## 4 healthy vs. no rejection D_365 6.321421e-09 1.896426e-08
## 5 healthy vs. acute rejection D_0 1.340189e-10 7.035992e-10
## 6 healthy vs. acute rejection D_7 2.975331e-10 1.249639e-09
## 7 healthy vs. acute rejection D_90 2.122261e-09 7.427913e-09
## 8 healthy vs. acute rejection D_365 9.383021e-06 2.189372e-05
## 9 healthy vs. acute rejection D_biopsy 1.355816e-07 3.559017e-07
## 10 no rejection D_0 vs. no rejection D_7 6.409870e-05 1.346073e-04
## 11 no rejection D_90 vs. no rejection D_365 2.581426e-02 4.517495e-02
## 12 acute rejection D_90 vs. acute rejection D_365 2.400159e-02 4.517495e-02
## star
## 1 ***
## 2 ***
## 3 ***
## 4 ***
## 5 ***
## 6 ***
## 7 ***
## 8 ***
## 9 ***
## 10 ***
## 11 *
## 12 *
pair_rest <- pair_table %>%
filter(p<0.05,
term1 != 'healthy.000') %>%
mutate(star = if_else(p<0.05, '*', '')) %>%
mutate(star = if_else(p<0.01, '**', star)) %>%
mutate(star = if_else(p<0.001, '***', star)) %>%
data.frame()3.3.3 Visualization
Open code
## plot
plotos <- data_long %>%
filter(!is.na(il1_ra_value)) %>%
mutate(group = factor(group, levels =
c("no_rejection",
"acute_rejection",
"healthy")),
time_point = factor(time_point)) %>%
mutate(time_point = if_else(group == 'healthy', 'healthy', time_point)) %>%
mutate(group_sep = interaction(group, time_point)) %>%
## plotting
ggplot(aes(x=time_point, y=il1_ra_value)) +
geom_boxplot(aes(fill=group, group = group_sep), outlier.shape = NA,
size=0.5,position = position_dodge(width = 0.84), color='black') +
labs(x = "Time post TX", y = "IL-1 RA (pg/ml)") +
scale_x_discrete(labels = c("Before TX", "7 days", "3 months", "365 days",
"Biopsy", "Healthy")) +
theme(axis.text = element_text(size = 11),
axis.title = element_text(size = 12),
legend.text = element_text(size = 11),
legend.title = element_text(size = 12)) +
scale_fill_manual(values = c("no_rejection" = "skyblue2",
"acute_rejection" = "coral",
"healthy" = "darkolivegreen3"),
labels = c("no rejection", "acute rejection", "healthy")) +
coord_cartesian(ylim = c(0, 7000)) +
annotate("text", x = 6, y = 2000, label = "#", size = 7)
pair_rest$y <- c(6900, 4000, 3600)
plotos <- plotos + geom_segment(data = pair_rest,
aes(x = xcoord_1, xend = xcoord_2, y = y, yend = y),
color = pair_rest$cole,
linetype = "solid", size = 0.4)
plot_il1_ra <- plotos + annotate("text", x = pair_rest$anot_coord, y = pair_rest$y+50,
label = pair_rest$star, size = 7, color=pair_rest$cole)
plot_il1_ra
3.4 il18
3.4.1 Mixed model
Open code
## filtering data
model_data <- data_long %>%
filter(!is.na(il18_value))
## model fit
model_il18 <- lmer(log10(il18_value) ~
group_sep +
(1|id), data = model_data)
## diagnostic plot of model
plot(model_il18)
Open code
## anova result od model
Anova(model_il18, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il18_value)
## F Df Df.res Pr(>F)
## group_sep 23.496 9 412.76 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## model summary
summary(model_il18)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il18_value) ~ group_sep + (1 | id)
## Data: model_data
##
## REML criterion at convergence: 290.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4547 -0.5610 -0.0786 0.4362 3.2320
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.03004 0.1733
## Residual 0.07579 0.2753
## Number of obs: 500, groups: id, 186
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.706285 0.051611 52.436
## group_sephealthy.000 -0.433108 0.069774 -6.207
## group_sepno_rejection.000 0.181152 0.061559 2.943
## group_sepacute_rejection.007 -0.250706 0.061125 -4.102
## group_sepno_rejection.007 -0.173250 0.061559 -2.814
## group_sepacute_rejection.090 0.133694 0.062099 2.153
## group_sepno_rejection.090 0.124231 0.071958 1.726
## group_sepacute_rejection.365 0.050040 0.065106 0.769
## group_sepno_rejection.365 -0.003025 0.071958 -0.042
## group_sepacute_rejection.biopsy 0.127269 0.069595 1.829
##
## Correlation of Fixed Effects:
## (Intr) g_.000 g__.000 grp_spc_.007 grp_spn_.007 grp_spc_.090
## grp_sph.000 -0.740
## grp_sp_.000 -0.838 0.620
## grp_spc_.007 -0.626 0.463 0.525
## grp_spn_.007 -0.838 0.620 0.787 0.525
## grp_spc_.090 -0.619 0.458 0.519 0.521 0.519
## grp_spn_.090 -0.717 0.531 0.673 0.449 0.673 0.444
## grp_spc_.365 -0.593 0.438 0.497 0.498 0.497 0.492
## grp_spn_.365 -0.717 0.531 0.673 0.449 0.673 0.444
## grp_spct_r. -0.551 0.407 0.462 0.462 0.462 0.458
## grp_spn_.090 grp_spc_.365 grp_spn_.365
## grp_sph.000
## grp_sp_.000
## grp_spc_.007
## grp_spn_.007
## grp_spc_.090
## grp_spn_.090
## grp_spc_.365 0.425
## grp_spn_.365 0.625 0.425
## grp_spct_r. 0.395 0.432 0.3953.4.2 Pairwise comparison
3.4.2.1 List of contrasts
List of contrast names and left and right coordinates on the X axis
Open code
contrast <- c(
## healthy vs. others
"healthy vs. no rejection D_0",
"healthy vs. no rejection D_7",
"healthy vs. no rejection D_90",
"healthy vs. no rejection D_365",
"healthy vs. acute rejection D_0",
"healthy vs. acute rejection D_7",
"healthy vs. acute rejection D_90",
"healthy vs. acute rejection D_365",
"healthy vs. acute rejection D_biopsy",
## between-subjects comparisons
"no rejection D_0 vs. acute rejection D_0",
"no rejection D_7 vs. acute rejection D_7",
"no rejection D_90 vs. acute rejection D_90",
"no rejection D_365 vs. acute rejection D_365",
"no rejection D_0 vs. acute rejection D_biopsy",
"no rejection D_7 vs. acute rejection D_biopsy",
"no rejection D_90 vs. acute rejection D_biopsy",
"no rejection D_365 vs. acute rejection D_biopsy",
## within-subjects comparisons
### no rejection
"no rejection D_0 vs. no rejection D_7",
"no rejection D_0 vs. no rejection D_90",
"no rejection D_0 vs. no rejection D_365",
"no rejection D_7 vs. no rejection D_90",
"no rejection D_7 vs. no rejection D_365",
"no rejection D_90 vs. no rejection D_365",
### acute rejection
"acute rejection D_0 vs. acute rejection D_7",
"acute rejection D_0 vs. acute rejection D_90",
"acute rejection D_0 vs. acute rejection D_365",
"acute rejection D_7 vs. acute rejection D_90",
"acute rejection D_7 vs. acute rejection D_365",
"acute rejection D_90 vs. acute rejection D_365",
"acute rejection D_0 vs. acute rejection D_biopsy",
"acute rejection D_7 vs. acute rejection D_biopsy",
"acute rejection D_90 vs. acute rejection D_biopsy",
"acute rejection D_365 vs. acute rejection D_biopsy"
)
xcoord_1 <- c(rep(6, 9), seq(0.8, 3.8, 1), seq(0.8, 3.8, 1),
rep(0.8, 3), rep(1.8, 2), 2.8,
rep(1.2, 3), rep(2.2, 2), 3.2,
seq(1.2, 4.2, 1))
xcoord_2 <- c(seq(0.8, 3.8, 1), seq(1.2, 4.2, 1),
5, seq(1.2, 4.2, 1), rep(5, 4),
seq(1.8, 3.8, 1), 2.8, 3.8, 3.8,
2.2, 3.2, 4.2, 3.2, 4.2, 4.2,
rep(5, 4))
anot_coord <- (xcoord_1+xcoord_2)/2
pair_table <- data.frame(contrast, xcoord_1, xcoord_2,
anot_coord)
pair_table <- pair_table %>%
mutate(term1 = if_else(xcoord_1 == 6, 'healthy.000', 'other')) %>%
mutate(term1 = if_else(xcoord_1 == 0.8, 'no_rejection.000', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 1.8, 'no_rejection.007', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 2.8, 'no_rejection.090', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 3.8, 'no_rejection.365', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 5, 'acute_rejection.biopsy', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 1.2, 'acute_rejection.000', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 2.2, 'acute_rejection.007', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 3.2, 'acute_rejection.090', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 4.2, 'acute_rejection.365', term1)) %>%
mutate(term2 = if_else(xcoord_2 == 6, 'healthy.000', 'other')) %>%
mutate(term2 = if_else(xcoord_2 == 0.8, 'no_rejection.000', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 1.8, 'no_rejection.007', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 2.8, 'no_rejection.090', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 3.8, 'no_rejection.365', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 5, 'acute_rejection.biopsy', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 1.2, 'acute_rejection.000', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 2.2, 'acute_rejection.007', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 3.2, 'acute_rejection.090', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 4.2, 'acute_rejection.365', term2))
pair_table$p <- NA3.4.2.2 Between-subjects comparisons
Open code
set.seed(16)
for (i in 1:17){
dat_pair <- model_data %>%
filter(group_sep == paste(pair_table$term1[i])|
group_sep == paste(pair_table$term2[i]))
g1 <- dat_pair %>% filter(group_sep == paste(pair_table$term1[i]))
g2 <- dat_pair %>% filter(group_sep == paste(pair_table$term2[i]))
pair_table$p[i] <- wilcox.test(g1$il18_value, g2$il18_value, paired=FALSE)$p.value
}3.4.2.3 Within-subject comparisons
Open code
set.seed(16)
for (i in 18:33){
dat_pair <- model_data %>%
filter(group_sep == paste(pair_table$term1[i]) |
group_sep == paste(pair_table$term2[i])) %>%
select(id, group_sep, il18_value) %>%
pivot_wider(names_from = group_sep, values_from = il18_value) %>%
filter(complete.cases(.)) %>% data.frame()
g1 <- dat_pair[,2]
g2 <- dat_pair[,3]
pair_table$p[i] <- wilcox.test(c(g1), c(g2), paired=TRUE)$p.value
}
pair_table <- pair_table[-c(14:15, 19:20, 22, 25:26, 28, 30:33),]
pair_table$cole <- c(rep("black", 15), rep("blue3", 3), rep("red4", 3))
pair_table$fdr <- p.adjust(pair_table$p, method = 'BH')3.4.2.4 Resulting table
Open code
## the whole `pair table`
pair_table
## contrast xcoord_1 xcoord_2 anot_coord
## 1 healthy vs. no rejection D_0 6.0 0.8 3.4
## 2 healthy vs. no rejection D_7 6.0 1.8 3.9
## 3 healthy vs. no rejection D_90 6.0 2.8 4.4
## 4 healthy vs. no rejection D_365 6.0 3.8 4.9
## 5 healthy vs. acute rejection D_0 6.0 1.2 3.6
## 6 healthy vs. acute rejection D_7 6.0 2.2 4.1
## 7 healthy vs. acute rejection D_90 6.0 3.2 4.6
## 8 healthy vs. acute rejection D_365 6.0 4.2 5.1
## 9 healthy vs. acute rejection D_biopsy 6.0 5.0 5.5
## 10 no rejection D_0 vs. acute rejection D_0 0.8 1.2 1.0
## 11 no rejection D_7 vs. acute rejection D_7 1.8 2.2 2.0
## 12 no rejection D_90 vs. acute rejection D_90 2.8 3.2 3.0
## 13 no rejection D_365 vs. acute rejection D_365 3.8 4.2 4.0
## 16 no rejection D_90 vs. acute rejection D_biopsy 2.8 5.0 3.9
## 17 no rejection D_365 vs. acute rejection D_biopsy 3.8 5.0 4.4
## 18 no rejection D_0 vs. no rejection D_7 0.8 1.8 1.3
## 21 no rejection D_7 vs. no rejection D_90 1.8 2.8 2.3
## 23 no rejection D_90 vs. no rejection D_365 2.8 3.8 3.3
## 24 acute rejection D_0 vs. acute rejection D_7 1.2 2.2 1.7
## 27 acute rejection D_7 vs. acute rejection D_90 2.2 3.2 2.7
## 29 acute rejection D_90 vs. acute rejection D_365 3.2 4.2 3.7
## term1 term2 p cole fdr
## 1 healthy.000 no_rejection.000 1.671896e-16 black 3.510982e-15
## 2 healthy.000 no_rejection.007 8.598524e-08 black 2.006322e-07
## 3 healthy.000 no_rejection.090 2.667314e-13 black 1.867120e-12
## 4 healthy.000 no_rejection.365 2.761390e-10 black 8.284170e-10
## 5 healthy.000 acute_rejection.000 2.124780e-06 black 4.056397e-06
## 6 healthy.000 acute_rejection.007 4.471962e-02 black 5.524189e-02
## 7 healthy.000 acute_rejection.090 8.245243e-13 black 4.328752e-12
## 8 healthy.000 acute_rejection.365 1.675587e-12 black 7.037463e-12
## 9 healthy.000 acute_rejection.biopsy 4.551076e-11 black 1.592877e-10
## 10 no_rejection.000 acute_rejection.000 2.664446e-02 black 3.497086e-02
## 11 no_rejection.007 acute_rejection.007 1.305168e-01 black 1.522696e-01
## 12 no_rejection.090 acute_rejection.090 7.088759e-01 black 7.088759e-01
## 13 no_rejection.365 acute_rejection.365 1.722459e-02 black 2.411443e-02
## 16 no_rejection.090 acute_rejection.biopsy 3.083806e-01 black 3.237997e-01
## 17 no_rejection.365 acute_rejection.biopsy 2.268561e-03 black 3.664598e-03
## 18 no_rejection.000 no_rejection.007 1.006443e-13 blue3 1.056765e-12
## 21 no_rejection.007 no_rejection.090 1.011576e-07 blue3 2.124310e-07
## 23 no_rejection.090 no_rejection.365 6.539723e-03 blue3 9.809584e-03
## 24 acute_rejection.000 acute_rejection.007 4.081365e-05 red4 7.142390e-05
## 27 acute_rejection.007 acute_rejection.090 3.163223e-09 red4 8.303459e-09
## 29 acute_rejection.090 acute_rejection.365 2.312297e-01 red4 2.555696e-01
## only significant differences
pair_table %>%
filter(p<0.05) %>%
select(contrast, p, fdr) %>%
mutate(star = if_else(p<0.05, '*', '')) %>%
mutate(star = if_else(p<0.01, '**', star)) %>%
mutate(star = if_else(p<0.001, '***', star)) %>%
data.frame()
## contrast p fdr
## 1 healthy vs. no rejection D_0 1.671896e-16 3.510982e-15
## 2 healthy vs. no rejection D_7 8.598524e-08 2.006322e-07
## 3 healthy vs. no rejection D_90 2.667314e-13 1.867120e-12
## 4 healthy vs. no rejection D_365 2.761390e-10 8.284170e-10
## 5 healthy vs. acute rejection D_0 2.124780e-06 4.056397e-06
## 6 healthy vs. acute rejection D_7 4.471962e-02 5.524189e-02
## 7 healthy vs. acute rejection D_90 8.245243e-13 4.328752e-12
## 8 healthy vs. acute rejection D_365 1.675587e-12 7.037463e-12
## 9 healthy vs. acute rejection D_biopsy 4.551076e-11 1.592877e-10
## 10 no rejection D_0 vs. acute rejection D_0 2.664446e-02 3.497086e-02
## 11 no rejection D_365 vs. acute rejection D_365 1.722459e-02 2.411443e-02
## 12 no rejection D_365 vs. acute rejection D_biopsy 2.268561e-03 3.664598e-03
## 13 no rejection D_0 vs. no rejection D_7 1.006443e-13 1.056765e-12
## 14 no rejection D_7 vs. no rejection D_90 1.011576e-07 2.124310e-07
## 15 no rejection D_90 vs. no rejection D_365 6.539723e-03 9.809584e-03
## 16 acute rejection D_0 vs. acute rejection D_7 4.081365e-05 7.142390e-05
## 17 acute rejection D_7 vs. acute rejection D_90 3.163223e-09 8.303459e-09
## star
## 1 ***
## 2 ***
## 3 ***
## 4 ***
## 5 ***
## 6 *
## 7 ***
## 8 ***
## 9 ***
## 10 *
## 11 *
## 12 **
## 13 ***
## 14 ***
## 15 **
## 16 ***
## 17 ***
pair_rest <- pair_table %>%
filter(p<0.05,
term1 != 'healthy.000') %>%
mutate(star = if_else(p<0.05, '*', '')) %>%
mutate(star = if_else(p<0.01, '**', star)) %>%
mutate(star = if_else(p<0.001, '***', star)) %>%
data.frame()3.4.3 Visualization
Open code
## plot
plotos <- data_long %>%
filter(!is.na(il18_value)) %>%
mutate(group = factor(group, levels =
c("no_rejection",
"acute_rejection",
"healthy")),
time_point = factor(time_point)) %>%
mutate(time_point = if_else(group == 'healthy', 'healthy', time_point)) %>%
mutate(group_sep = interaction(group, time_point)) %>%
## plotting
ggplot(aes(x=time_point, y=il18_value)) +
geom_boxplot(aes(fill=group, group = group_sep), outlier.shape = NA,
size=0.5,position = position_dodge(width = 0.84), color='black') +
labs(x = "Time post TX", y = "IL-18 (pg/ml)") +
scale_x_discrete(labels = c("Before TX", "7 days", "3 months", "365 days",
"Biopsy", "Healthy")) +
theme(axis.text = element_text(size = 11),
axis.title = element_text(size = 12),
legend.text = element_text(size = 11),
legend.title = element_text(size = 12)) +
scale_fill_manual(values = c("no_rejection" = "skyblue2",
"acute_rejection" = "coral",
"healthy" = "darkolivegreen3"),
labels = c("no rejection", "acute rejection", "healthy")) +
coord_cartesian(ylim = c(0, 3000)) +
annotate("text", x = 6, y = 700, label = "#", size = 7)
pair_rest$y <- c(2700, 1000, 1700, 2900, 2300, 2100, 1900, 1700)
plotos <- plotos + geom_segment(data = pair_rest,
aes(x = xcoord_1, xend = xcoord_2, y = y, yend = y),
color = pair_rest$cole,
linetype = "solid", size = 0.4)
plot_il18 <- plotos + annotate("text", x = pair_rest$anot_coord, y = pair_rest$y+50,
label = pair_rest$star, size = 7, color=pair_rest$cole)
plot_il18
3.5 il18_bp
3.5.1 Mixed model
Open code
## filtering data
model_data <- data_long %>%
filter(!is.na(il18_bp_value))
## model fit
model_il18_bp <- lmer(log10(il18_bp_value) ~
group_sep +
(1|id), data = model_data)
## diagnostic plot of model
plot(model_il18_bp)
Open code
## anova result od model
Anova(model_il18_bp, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il18_bp_value)
## F Df Df.res Pr(>F)
## group_sep 8.4695 9 489.8 8.293e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## model summary
summary(model_il18_bp)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il18_bp_value) ~ group_sep + (1 | id)
## Data: model_data
##
## REML criterion at convergence: -484.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -10.0578 -0.3070 0.0384 0.3462 4.8019
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.006046 0.07775
## Residual 0.018843 0.13727
## Number of obs: 577, groups: id, 169
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.614819 0.024804 145.735
## group_sephealthy.000 -0.165875 0.037658 -4.405
## group_sepno_rejection.000 0.004244 0.029707 0.143
## group_sepacute_rejection.007 -0.061525 0.030101 -2.044
## group_sepno_rejection.007 -0.046670 0.029665 -1.573
## group_sepacute_rejection.090 -0.121220 0.030711 -3.947
## group_sepno_rejection.090 -0.109015 0.029982 -3.636
## group_sepacute_rejection.365 -0.092007 0.031914 -2.883
## group_sepno_rejection.365 -0.121527 0.030184 -4.026
## group_sepacute_rejection.biopsy -0.053131 0.034328 -1.548
##
## Correlation of Fixed Effects:
## (Intr) g_.000 g__.000 grp_spc_.007 grp_spn_.007 grp_spc_.090
## grp_sph.000 -0.659
## grp_sp_.000 -0.835 0.550
## grp_spc_.007 -0.640 0.422 0.534
## grp_spn_.007 -0.836 0.551 0.771 0.535
## grp_spc_.090 -0.628 0.414 0.524 0.517 0.525
## grp_spn_.090 -0.827 0.545 0.763 0.529 0.764 0.520
## grp_spc_.365 -0.606 0.399 0.506 0.499 0.507 0.489
## grp_spn_.365 -0.822 0.541 0.758 0.526 0.759 0.516
## grp_spct_r. -0.556 0.366 0.464 0.458 0.465 0.449
## grp_spn_.090 grp_spc_.365 grp_spn_.365
## grp_sph.000
## grp_sp_.000
## grp_spc_.007
## grp_spn_.007
## grp_spc_.090
## grp_spn_.090
## grp_spc_.365 0.501
## grp_spn_.365 0.752 0.498
## grp_spct_r. 0.460 0.428 0.4573.5.2 Pairwise comparison
3.5.2.1 List of contrasts
List of contrast names and left and right coordinates on the X axis
Open code
contrast <- c(
## healthy vs. others
"healthy vs. no rejection D_0",
"healthy vs. no rejection D_7",
"healthy vs. no rejection D_90",
"healthy vs. no rejection D_365",
"healthy vs. acute rejection D_0",
"healthy vs. acute rejection D_7",
"healthy vs. acute rejection D_90",
"healthy vs. acute rejection D_365",
"healthy vs. acute rejection D_biopsy",
## between-subjects comparisons
"no rejection D_0 vs. acute rejection D_0",
"no rejection D_7 vs. acute rejection D_7",
"no rejection D_90 vs. acute rejection D_90",
"no rejection D_365 vs. acute rejection D_365",
"no rejection D_0 vs. acute rejection D_biopsy",
"no rejection D_7 vs. acute rejection D_biopsy",
"no rejection D_90 vs. acute rejection D_biopsy",
"no rejection D_365 vs. acute rejection D_biopsy",
## within-subjects comparisons
### no rejection
"no rejection D_0 vs. no rejection D_7",
"no rejection D_0 vs. no rejection D_90",
"no rejection D_0 vs. no rejection D_365",
"no rejection D_7 vs. no rejection D_90",
"no rejection D_7 vs. no rejection D_365",
"no rejection D_90 vs. no rejection D_365",
### acute rejection
"acute rejection D_0 vs. acute rejection D_7",
"acute rejection D_0 vs. acute rejection D_90",
"acute rejection D_0 vs. acute rejection D_365",
"acute rejection D_7 vs. acute rejection D_90",
"acute rejection D_7 vs. acute rejection D_365",
"acute rejection D_90 vs. acute rejection D_365",
"acute rejection D_0 vs. acute rejection D_biopsy",
"acute rejection D_7 vs. acute rejection D_biopsy",
"acute rejection D_90 vs. acute rejection D_biopsy",
"acute rejection D_365 vs. acute rejection D_biopsy"
)
xcoord_1 <- c(rep(6, 9), seq(0.8, 3.8, 1), seq(0.8, 3.8, 1),
rep(0.8, 3), rep(1.8, 2), 2.8,
rep(1.2, 3), rep(2.2, 2), 3.2,
seq(1.2, 4.2, 1))
xcoord_2 <- c(seq(0.8, 3.8, 1), seq(1.2, 4.2, 1),
5, seq(1.2, 4.2, 1), rep(5, 4),
seq(1.8, 3.8, 1), 2.8, 3.8, 3.8,
2.2, 3.2, 4.2, 3.2, 4.2, 4.2,
rep(5, 4))
anot_coord <- (xcoord_1+xcoord_2)/2
pair_table <- data.frame(contrast, xcoord_1, xcoord_2,
anot_coord)
pair_table <- pair_table %>%
mutate(term1 = if_else(xcoord_1 == 6, 'healthy.000', 'other')) %>%
mutate(term1 = if_else(xcoord_1 == 0.8, 'no_rejection.000', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 1.8, 'no_rejection.007', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 2.8, 'no_rejection.090', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 3.8, 'no_rejection.365', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 5, 'acute_rejection.biopsy', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 1.2, 'acute_rejection.000', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 2.2, 'acute_rejection.007', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 3.2, 'acute_rejection.090', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 4.2, 'acute_rejection.365', term1)) %>%
mutate(term2 = if_else(xcoord_2 == 6, 'healthy.000', 'other')) %>%
mutate(term2 = if_else(xcoord_2 == 0.8, 'no_rejection.000', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 1.8, 'no_rejection.007', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 2.8, 'no_rejection.090', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 3.8, 'no_rejection.365', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 5, 'acute_rejection.biopsy', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 1.2, 'acute_rejection.000', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 2.2, 'acute_rejection.007', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 3.2, 'acute_rejection.090', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 4.2, 'acute_rejection.365', term2))
pair_table$p <- NA3.5.2.2 Between-subjects comparisons
Open code
set.seed(16)
for (i in 1:17){
dat_pair <- model_data %>%
filter(group_sep == paste(pair_table$term1[i])|
group_sep == paste(pair_table$term2[i]))
g1 <- dat_pair %>% filter(group_sep == paste(pair_table$term1[i]))
g2 <- dat_pair %>% filter(group_sep == paste(pair_table$term2[i]))
pair_table$p[i] <- wilcox.test(g1$il18_bp_value, g2$il18_bp_value, paired=FALSE)$p.value
}3.5.2.3 Within-subject comparisons
Open code
set.seed(16)
for (i in 18:33){
dat_pair <- model_data %>%
filter(group_sep == paste(pair_table$term1[i]) |
group_sep == paste(pair_table$term2[i])) %>%
select(id, group_sep, il18_bp_value) %>%
pivot_wider(names_from = group_sep, values_from = il18_bp_value) %>%
filter(complete.cases(.)) %>% data.frame()
g1 <- dat_pair[,2]
g2 <- dat_pair[,3]
pair_table$p[i] <- wilcox.test(c(g1), c(g2), paired=TRUE)$p.value
}
pair_table <- pair_table[-c(14:15, 19:20, 22, 25:26, 28, 30:33),]
pair_table$cole <- c(rep("black", 15), rep("blue3", 3), rep("red4", 3))
pair_table$fdr <- p.adjust(pair_table$p, method = 'BH')3.5.2.4 Resulting table
Open code
## the whole `pair table`
pair_table
## contrast xcoord_1 xcoord_2 anot_coord
## 1 healthy vs. no rejection D_0 6.0 0.8 3.4
## 2 healthy vs. no rejection D_7 6.0 1.8 3.9
## 3 healthy vs. no rejection D_90 6.0 2.8 4.4
## 4 healthy vs. no rejection D_365 6.0 3.8 4.9
## 5 healthy vs. acute rejection D_0 6.0 1.2 3.6
## 6 healthy vs. acute rejection D_7 6.0 2.2 4.1
## 7 healthy vs. acute rejection D_90 6.0 3.2 4.6
## 8 healthy vs. acute rejection D_365 6.0 4.2 5.1
## 9 healthy vs. acute rejection D_biopsy 6.0 5.0 5.5
## 10 no rejection D_0 vs. acute rejection D_0 0.8 1.2 1.0
## 11 no rejection D_7 vs. acute rejection D_7 1.8 2.2 2.0
## 12 no rejection D_90 vs. acute rejection D_90 2.8 3.2 3.0
## 13 no rejection D_365 vs. acute rejection D_365 3.8 4.2 4.0
## 16 no rejection D_90 vs. acute rejection D_biopsy 2.8 5.0 3.9
## 17 no rejection D_365 vs. acute rejection D_biopsy 3.8 5.0 4.4
## 18 no rejection D_0 vs. no rejection D_7 0.8 1.8 1.3
## 21 no rejection D_7 vs. no rejection D_90 1.8 2.8 2.3
## 23 no rejection D_90 vs. no rejection D_365 2.8 3.8 3.3
## 24 acute rejection D_0 vs. acute rejection D_7 1.2 2.2 1.7
## 27 acute rejection D_7 vs. acute rejection D_90 2.2 3.2 2.7
## 29 acute rejection D_90 vs. acute rejection D_365 3.2 4.2 3.7
## term1 term2 p cole fdr
## 1 healthy.000 no_rejection.000 4.111120e-12 black 8.633353e-11
## 2 healthy.000 no_rejection.007 1.240228e-08 black 8.681598e-08
## 3 healthy.000 no_rejection.090 2.278600e-03 black 4.350054e-03
## 4 healthy.000 no_rejection.365 1.381490e-04 black 3.626412e-04
## 5 healthy.000 acute_rejection.000 1.711485e-09 black 1.797059e-08
## 6 healthy.000 acute_rejection.007 9.191676e-05 black 2.757503e-04
## 7 healthy.000 acute_rejection.090 4.719442e-02 black 6.607219e-02
## 8 healthy.000 acute_rejection.365 8.426347e-03 black 1.361179e-02
## 9 healthy.000 acute_rejection.biopsy 1.333362e-06 black 5.600122e-06
## 10 no_rejection.000 acute_rejection.000 3.168633e-01 black 3.696738e-01
## 11 no_rejection.007 acute_rejection.007 1.275900e-01 black 1.674619e-01
## 12 no_rejection.090 acute_rejection.090 4.543678e-01 black 5.021959e-01
## 13 no_rejection.365 acute_rejection.365 7.350379e-01 black 7.350379e-01
## 16 no_rejection.090 acute_rejection.biopsy 4.632816e-03 black 8.107427e-03
## 17 no_rejection.365 acute_rejection.biopsy 3.163257e-02 black 4.744886e-02
## 18 no_rejection.000 no_rejection.007 2.090599e-05 blue3 7.317096e-05
## 21 no_rejection.007 no_rejection.090 4.544604e-08 blue3 2.385917e-07
## 23 no_rejection.090 no_rejection.365 5.092835e-01 blue3 5.347477e-01
## 24 acute_rejection.000 acute_rejection.007 1.367483e-03 red4 2.943894e-03
## 27 acute_rejection.007 acute_rejection.090 1.401854e-03 red4 2.943894e-03
## 29 acute_rejection.090 acute_rejection.365 1.547263e-01 red4 1.911325e-01
## only significant differences
pair_table %>%
filter(p<0.05) %>%
select(contrast, p, fdr) %>%
mutate(star = if_else(p<0.05, '*', '')) %>%
mutate(star = if_else(p<0.01, '**', star)) %>%
mutate(star = if_else(p<0.001, '***', star)) %>%
data.frame()
## contrast p fdr
## 1 healthy vs. no rejection D_0 4.111120e-12 8.633353e-11
## 2 healthy vs. no rejection D_7 1.240228e-08 8.681598e-08
## 3 healthy vs. no rejection D_90 2.278600e-03 4.350054e-03
## 4 healthy vs. no rejection D_365 1.381490e-04 3.626412e-04
## 5 healthy vs. acute rejection D_0 1.711485e-09 1.797059e-08
## 6 healthy vs. acute rejection D_7 9.191676e-05 2.757503e-04
## 7 healthy vs. acute rejection D_90 4.719442e-02 6.607219e-02
## 8 healthy vs. acute rejection D_365 8.426347e-03 1.361179e-02
## 9 healthy vs. acute rejection D_biopsy 1.333362e-06 5.600122e-06
## 10 no rejection D_90 vs. acute rejection D_biopsy 4.632816e-03 8.107427e-03
## 11 no rejection D_365 vs. acute rejection D_biopsy 3.163257e-02 4.744886e-02
## 12 no rejection D_0 vs. no rejection D_7 2.090599e-05 7.317096e-05
## 13 no rejection D_7 vs. no rejection D_90 4.544604e-08 2.385917e-07
## 14 acute rejection D_0 vs. acute rejection D_7 1.367483e-03 2.943894e-03
## 15 acute rejection D_7 vs. acute rejection D_90 1.401854e-03 2.943894e-03
## star
## 1 ***
## 2 ***
## 3 **
## 4 ***
## 5 ***
## 6 ***
## 7 *
## 8 **
## 9 ***
## 10 **
## 11 *
## 12 ***
## 13 ***
## 14 **
## 15 **
pair_rest <- pair_table %>%
filter(p<0.05,
term1 != 'healthy.000') %>%
mutate(star = if_else(p<0.05, '*', '')) %>%
mutate(star = if_else(p<0.01, '**', star)) %>%
mutate(star = if_else(p<0.001, '***', star)) %>%
data.frame()! Healthy group does not differ from others with p<0.001, but some only p<0.05
3.5.3 Visualization
Open code
## plot
plotos <- data_long %>%
filter(!is.na(il18_bp_value)) %>%
mutate(group = factor(group, levels =
c("no_rejection",
"acute_rejection",
"healthy")),
time_point = factor(time_point)) %>%
mutate(time_point = if_else(group == 'healthy', 'healthy', time_point)) %>%
mutate(group_sep = interaction(group, time_point)) %>%
## plotting
ggplot(aes(x=time_point, y=il18_bp_value)) +
geom_boxplot(aes(fill=group, group = group_sep), outlier.shape = NA,
size=0.5,position = position_dodge(width = 0.84), color='black') +
labs(x = "Time post TX", y = "IL-18 BP (pg/ml)") +
scale_x_discrete(labels = c("Before TX", "7 days", "3 months", "365 days",
"Biopsy", "Healthy")) +
theme(axis.text = element_text(size = 11),
axis.title = element_text(size = 12),
legend.text = element_text(size = 11),
legend.title = element_text(size = 12)) +
scale_fill_manual(values = c("no_rejection" = "skyblue2",
"acute_rejection" = "coral",
"healthy" = "darkolivegreen3"),
labels = c("no rejection", "acute rejection", "healthy")) +
coord_cartesian(ylim = c(0, 7300)) +
annotate("text", x = 6, y = 5000, label = "#", size = 7)
pair_rest$y <- c(5800, 5500, 7200, 6900, 6400, 6100)
plotos <- plotos + geom_segment(data = pair_rest,
aes(x = xcoord_1, xend = xcoord_2, y = y, yend = y),
color = pair_rest$cole,
linetype = "solid", size = 0.4)
plot_il18_bp <- plotos + annotate("text", x = pair_rest$anot_coord, y = pair_rest$y+50,
label = pair_rest$star, size = 7, color=pair_rest$cole)
plot_il18_bp
3.6 il18_free
3.6.1 Mixed model
Open code
## filtering data
model_data <- data_long %>%
filter(!is.na(il18_free))
## model fit
model_il18_free <- lmer(log10(il18_free) ~
group_sep +
(1|id), data = model_data)
## diagnostic plot of model
plot(model_il18_free)
Open code
## anova result od model
Anova(model_il18_free, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il18_free)
## F Df Df.res Pr(>F)
## group_sep 19.077 9 401.57 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## model summary
summary(model_il18_free)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il18_free) ~ group_sep + (1 | id)
## Data: model_data
##
## REML criterion at convergence: 247.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4008 -0.5625 -0.0955 0.4128 3.3292
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.02494 0.1579
## Residual 0.07261 0.2695
## Number of obs: 481, groups: id, 169
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.51490 0.05020 50.101
## group_sephealthy.000 -0.33233 0.07528 -4.415
## group_sepno_rejection.000 0.16742 0.05973 2.803
## group_sepacute_rejection.007 -0.23114 0.06029 -3.834
## group_sepno_rejection.007 -0.16130 0.05964 -2.704
## group_sepacute_rejection.090 0.17291 0.06123 2.824
## group_sepno_rejection.090 0.14789 0.06978 2.119
## group_sepacute_rejection.365 0.07912 0.06415 1.233
## group_sepno_rejection.365 0.01916 0.06978 0.275
## group_sepacute_rejection.biopsy 0.14763 0.06836 2.160
##
## Correlation of Fixed Effects:
## (Intr) g_.000 g__.000 grp_spc_.007 grp_spn_.007 grp_spc_.090
## grp_sph.000 -0.667
## grp_sp_.000 -0.840 0.560
## grp_spc_.007 -0.643 0.429 0.540
## grp_spn_.007 -0.842 0.561 0.782 0.541
## grp_spc_.090 -0.636 0.424 0.535 0.528 0.536
## grp_spn_.090 -0.719 0.480 0.669 0.463 0.669 0.458
## grp_spc_.365 -0.610 0.407 0.513 0.506 0.514 0.500
## grp_spn_.365 -0.719 0.480 0.669 0.463 0.669 0.458
## grp_spct_r. -0.566 0.378 0.476 0.469 0.477 0.465
## grp_spn_.090 grp_spc_.365 grp_spn_.365
## grp_sph.000
## grp_sp_.000
## grp_spc_.007
## grp_spn_.007
## grp_spc_.090
## grp_spn_.090
## grp_spc_.365 0.439
## grp_spn_.365 0.618 0.439
## grp_spct_r. 0.407 0.439 0.4073.6.2 Pairwise comparison
3.6.2.1 List of contrasts
List of contrast names and left and right coordinates on the X axis
Open code
contrast <- c(
## healthy vs. others
"healthy vs. no rejection D_0",
"healthy vs. no rejection D_7",
"healthy vs. no rejection D_90",
"healthy vs. no rejection D_365",
"healthy vs. acute rejection D_0",
"healthy vs. acute rejection D_7",
"healthy vs. acute rejection D_90",
"healthy vs. acute rejection D_365",
"healthy vs. acute rejection D_biopsy",
## between-subjects comparisons
"no rejection D_0 vs. acute rejection D_0",
"no rejection D_7 vs. acute rejection D_7",
"no rejection D_90 vs. acute rejection D_90",
"no rejection D_365 vs. acute rejection D_365",
"no rejection D_0 vs. acute rejection D_biopsy",
"no rejection D_7 vs. acute rejection D_biopsy",
"no rejection D_90 vs. acute rejection D_biopsy",
"no rejection D_365 vs. acute rejection D_biopsy",
## within-subjects comparisons
### no rejection
"no rejection D_0 vs. no rejection D_7",
"no rejection D_0 vs. no rejection D_90",
"no rejection D_0 vs. no rejection D_365",
"no rejection D_7 vs. no rejection D_90",
"no rejection D_7 vs. no rejection D_365",
"no rejection D_90 vs. no rejection D_365",
### acute rejection
"acute rejection D_0 vs. acute rejection D_7",
"acute rejection D_0 vs. acute rejection D_90",
"acute rejection D_0 vs. acute rejection D_365",
"acute rejection D_7 vs. acute rejection D_90",
"acute rejection D_7 vs. acute rejection D_365",
"acute rejection D_90 vs. acute rejection D_365",
"acute rejection D_0 vs. acute rejection D_biopsy",
"acute rejection D_7 vs. acute rejection D_biopsy",
"acute rejection D_90 vs. acute rejection D_biopsy",
"acute rejection D_365 vs. acute rejection D_biopsy"
)
xcoord_1 <- c(rep(6, 9), seq(0.8, 3.8, 1), seq(0.8, 3.8, 1),
rep(0.8, 3), rep(1.8, 2), 2.8,
rep(1.2, 3), rep(2.2, 2), 3.2,
seq(1.2, 4.2, 1))
xcoord_2 <- c(seq(0.8, 3.8, 1), seq(1.2, 4.2, 1),
5, seq(1.2, 4.2, 1), rep(5, 4),
seq(1.8, 3.8, 1), 2.8, 3.8, 3.8,
2.2, 3.2, 4.2, 3.2, 4.2, 4.2,
rep(5, 4))
anot_coord <- (xcoord_1+xcoord_2)/2
pair_table <- data.frame(contrast, xcoord_1, xcoord_2,
anot_coord)
pair_table <- pair_table %>%
mutate(term1 = if_else(xcoord_1 == 6, 'healthy.000', 'other')) %>%
mutate(term1 = if_else(xcoord_1 == 0.8, 'no_rejection.000', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 1.8, 'no_rejection.007', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 2.8, 'no_rejection.090', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 3.8, 'no_rejection.365', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 5, 'acute_rejection.biopsy', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 1.2, 'acute_rejection.000', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 2.2, 'acute_rejection.007', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 3.2, 'acute_rejection.090', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 4.2, 'acute_rejection.365', term1)) %>%
mutate(term2 = if_else(xcoord_2 == 6, 'healthy.000', 'other')) %>%
mutate(term2 = if_else(xcoord_2 == 0.8, 'no_rejection.000', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 1.8, 'no_rejection.007', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 2.8, 'no_rejection.090', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 3.8, 'no_rejection.365', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 5, 'acute_rejection.biopsy', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 1.2, 'acute_rejection.000', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 2.2, 'acute_rejection.007', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 3.2, 'acute_rejection.090', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 4.2, 'acute_rejection.365', term2))
pair_table$p <- NA3.6.2.2 Between-subjects comparisons
Open code
set.seed(16)
for (i in 1:17){
dat_pair <- model_data %>%
filter(group_sep == paste(pair_table$term1[i])|
group_sep == paste(pair_table$term2[i]))
g1 <- dat_pair %>% filter(group_sep == paste(pair_table$term1[i]))
g2 <- dat_pair %>% filter(group_sep == paste(pair_table$term2[i]))
pair_table$p[i] <- wilcox.test(g1$il18_free, g2$il18_free, paired=FALSE)$p.value
}3.6.2.3 Within-subject comparisons
Open code
set.seed(16)
for (i in 18:33){
dat_pair <- model_data %>%
filter(group_sep == paste(pair_table$term1[i]) |
group_sep == paste(pair_table$term2[i])) %>%
select(id, group_sep, il18_free) %>%
pivot_wider(names_from = group_sep, values_from = il18_free) %>%
filter(complete.cases(.)) %>% data.frame()
g1 <- dat_pair[,2]
g2 <- dat_pair[,3]
pair_table$p[i] <- wilcox.test(c(g1), c(g2), paired=TRUE)$p.value
}
pair_table <- pair_table[-c(14:15, 19:20, 22, 25:26, 28, 30:33),]
pair_table$cole <- c(rep("black", 15), rep("blue3", 3), rep("red4", 3))
pair_table$fdr <- p.adjust(pair_table$p, method = 'BH')3.6.2.4 Resulting table
Open code
## the whole `pair table`
pair_table
## contrast xcoord_1 xcoord_2 anot_coord
## 1 healthy vs. no rejection D_0 6.0 0.8 3.4
## 2 healthy vs. no rejection D_7 6.0 1.8 3.9
## 3 healthy vs. no rejection D_90 6.0 2.8 4.4
## 4 healthy vs. no rejection D_365 6.0 3.8 4.9
## 5 healthy vs. acute rejection D_0 6.0 1.2 3.6
## 6 healthy vs. acute rejection D_7 6.0 2.2 4.1
## 7 healthy vs. acute rejection D_90 6.0 3.2 4.6
## 8 healthy vs. acute rejection D_365 6.0 4.2 5.1
## 9 healthy vs. acute rejection D_biopsy 6.0 5.0 5.5
## 10 no rejection D_0 vs. acute rejection D_0 0.8 1.2 1.0
## 11 no rejection D_7 vs. acute rejection D_7 1.8 2.2 2.0
## 12 no rejection D_90 vs. acute rejection D_90 2.8 3.2 3.0
## 13 no rejection D_365 vs. acute rejection D_365 3.8 4.2 4.0
## 16 no rejection D_90 vs. acute rejection D_biopsy 2.8 5.0 3.9
## 17 no rejection D_365 vs. acute rejection D_biopsy 3.8 5.0 4.4
## 18 no rejection D_0 vs. no rejection D_7 0.8 1.8 1.3
## 21 no rejection D_7 vs. no rejection D_90 1.8 2.8 2.3
## 23 no rejection D_90 vs. no rejection D_365 2.8 3.8 3.3
## 24 acute rejection D_0 vs. acute rejection D_7 1.2 2.2 1.7
## 27 acute rejection D_7 vs. acute rejection D_90 2.2 3.2 2.7
## 29 acute rejection D_90 vs. acute rejection D_365 3.2 4.2 3.7
## term1 term2 p cole fdr
## 1 healthy.000 no_rejection.000 1.195713e-10 black 4.184996e-10
## 2 healthy.000 no_rejection.007 1.125705e-04 black 2.363980e-04
## 3 healthy.000 no_rejection.090 1.142557e-13 black 2.382870e-12
## 4 healthy.000 no_rejection.365 2.897492e-09 black 7.605918e-09
## 5 healthy.000 acute_rejection.000 1.390081e-04 black 2.646890e-04
## 6 healthy.000 acute_rejection.007 2.848719e-01 black 3.148584e-01
## 7 healthy.000 acute_rejection.090 2.269400e-13 black 2.382870e-12
## 8 healthy.000 acute_rejection.365 1.488787e-12 black 1.042151e-11
## 9 healthy.000 acute_rejection.biopsy 1.051022e-10 black 4.184996e-10
## 10 no_rejection.000 acute_rejection.000 3.508526e-02 black 4.604940e-02
## 11 no_rejection.007 acute_rejection.007 9.916951e-02 black 1.225035e-01
## 12 no_rejection.090 acute_rejection.090 4.724356e-01 black 4.724356e-01
## 13 no_rejection.365 acute_rejection.365 1.561394e-02 black 2.185951e-02
## 16 no_rejection.090 acute_rejection.biopsy 3.794611e-01 black 3.984342e-01
## 17 no_rejection.365 acute_rejection.biopsy 3.232913e-03 black 5.222398e-03
## 18 no_rejection.000 no_rejection.007 6.418077e-12 blue3 3.369490e-11
## 21 no_rejection.007 no_rejection.090 3.163586e-08 blue3 7.381702e-08
## 23 no_rejection.090 no_rejection.365 5.187312e-03 blue3 7.780968e-03
## 24 acute_rejection.000 acute_rejection.007 1.512508e-04 red4 2.646890e-04
## 27 acute_rejection.007 acute_rejection.090 1.644366e-09 red4 4.933099e-09
## 29 acute_rejection.090 acute_rejection.365 1.776244e-01 red4 2.072284e-01
## only significant differences
pair_table %>%
filter(p<0.05) %>%
select(contrast, p, fdr) %>%
mutate(star = if_else(p<0.05, '*', '')) %>%
mutate(star = if_else(p<0.01, '**', star)) %>%
mutate(star = if_else(p<0.001, '***', star)) %>%
data.frame()
## contrast p fdr
## 1 healthy vs. no rejection D_0 1.195713e-10 4.184996e-10
## 2 healthy vs. no rejection D_7 1.125705e-04 2.363980e-04
## 3 healthy vs. no rejection D_90 1.142557e-13 2.382870e-12
## 4 healthy vs. no rejection D_365 2.897492e-09 7.605918e-09
## 5 healthy vs. acute rejection D_0 1.390081e-04 2.646890e-04
## 6 healthy vs. acute rejection D_90 2.269400e-13 2.382870e-12
## 7 healthy vs. acute rejection D_365 1.488787e-12 1.042151e-11
## 8 healthy vs. acute rejection D_biopsy 1.051022e-10 4.184996e-10
## 9 no rejection D_0 vs. acute rejection D_0 3.508526e-02 4.604940e-02
## 10 no rejection D_365 vs. acute rejection D_365 1.561394e-02 2.185951e-02
## 11 no rejection D_365 vs. acute rejection D_biopsy 3.232913e-03 5.222398e-03
## 12 no rejection D_0 vs. no rejection D_7 6.418077e-12 3.369490e-11
## 13 no rejection D_7 vs. no rejection D_90 3.163586e-08 7.381702e-08
## 14 no rejection D_90 vs. no rejection D_365 5.187312e-03 7.780968e-03
## 15 acute rejection D_0 vs. acute rejection D_7 1.512508e-04 2.646890e-04
## 16 acute rejection D_7 vs. acute rejection D_90 1.644366e-09 4.933099e-09
## star
## 1 ***
## 2 ***
## 3 ***
## 4 ***
## 5 ***
## 6 ***
## 7 ***
## 8 ***
## 9 *
## 10 *
## 11 **
## 12 ***
## 13 ***
## 14 **
## 15 ***
## 16 ***
pair_rest <- pair_table %>%
filter(p<0.05,
term1 != 'healthy.000') %>%
mutate(star = if_else(p<0.05, '*', '')) %>%
mutate(star = if_else(p<0.01, '**', star)) %>%
mutate(star = if_else(p<0.001, '***', star)) %>%
data.frame()3.6.3 Visualization
Open code
## plot
plotos <- data_long %>%
filter(!is.na(il18_free)) %>%
mutate(group = factor(group, levels =
c("no_rejection",
"acute_rejection",
"healthy")),
time_point = factor(time_point)) %>%
mutate(time_point = if_else(group == 'healthy', 'healthy', time_point)) %>%
mutate(group_sep = interaction(group, time_point)) %>%
## plotting
ggplot(aes(x=time_point, y=il18_free)) +
geom_boxplot(aes(fill=group, group = group_sep), outlier.shape = NA,
size=0.5,position = position_dodge(width = 0.84), color='black') +
labs(x = "Time post TX", y = "Free IL-18 (pg/ml)") +
scale_x_discrete(labels = c("Before TX", "7 days", "3 months", "365 days",
"Biopsy", "Healthy")) +
theme(axis.text = element_text(size = 11),
axis.title = element_text(size = 12),
legend.text = element_text(size = 11),
legend.title = element_text(size = 12)) +
scale_fill_manual(values = c("no_rejection" = "skyblue2",
"acute_rejection" = "coral",
"healthy" = "darkolivegreen3"),
labels = c("no rejection", "acute rejection", "healthy")) +
coord_cartesian(ylim = c(0, 1700)) +
annotate("text", x = 6, y = 500, label = "#", size = 7)
pair_rest$y <- c(1650, 700, 1060, 1720,
830, 1200, 1350, 1280)
plotos <- plotos + geom_segment(data = pair_rest,
aes(x = xcoord_1, xend = xcoord_2, y = y, yend = y),
color = pair_rest$cole,
linetype = "solid", size = 0.4)
plot_il18_free <- plotos + annotate("text", x = pair_rest$anot_coord, y = pair_rest$y+5,
label = pair_rest$star, size = 7, color=pair_rest$cole)
plot_il18_free
3.7 il36_b
3.7.1 Mixed model
Open code
## filtering data
model_data <- data_long %>%
filter(!is.na(il36_b_value))
## model fit
model_il36_b <- lmer(log10(il36_b_value) ~
group_sep +
(1|id), data = model_data)
## diagnostic plot of model
plot(model_il36_b)
Open code
## anova result od model
Anova(model_il36_b, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il36_b_value)
## F Df Df.res Pr(>F)
## group_sep 13.599 9 411.76 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## model summary
summary(model_il36_b)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il36_b_value) ~ group_sep + (1 | id)
## Data: model_data
##
## REML criterion at convergence: 566.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.3627 -0.3397 0.0080 0.4167 3.2126
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.05778 0.2404
## Residual 0.13031 0.3610
## Number of obs: 500, groups: id, 186
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.70905 0.06872 10.318
## group_sephealthy.000 -0.52054 0.09296 -5.600
## group_sepno_rejection.000 0.05180 0.08200 0.632
## group_sepacute_rejection.007 -0.16912 0.08016 -2.110
## group_sepno_rejection.007 -0.05161 0.08200 -0.629
## group_sepacute_rejection.090 -0.36787 0.08146 -4.516
## group_sepno_rejection.090 -0.13766 0.09559 -1.440
## group_sepacute_rejection.365 -0.52969 0.08543 -6.200
## group_sepno_rejection.365 -0.19032 0.09559 -1.991
## group_sepacute_rejection.biopsy -0.50731 0.09134 -5.554
##
## Correlation of Fixed Effects:
## (Intr) g_.000 g__.000 grp_spc_.007 grp_spn_.007 grp_spc_.090
## grp_sph.000 -0.739
## grp_sp_.000 -0.838 0.620
## grp_spc_.007 -0.616 0.456 0.517
## grp_spn_.007 -0.838 0.620 0.794 0.517
## grp_spc_.090 -0.610 0.451 0.511 0.521 0.511
## grp_spn_.090 -0.719 0.531 0.681 0.443 0.681 0.439
## grp_spc_.365 -0.584 0.432 0.489 0.498 0.489 0.492
## grp_spn_.365 -0.719 0.531 0.681 0.443 0.681 0.439
## grp_spct_r. -0.542 0.401 0.455 0.462 0.455 0.458
## grp_spn_.090 grp_spc_.365 grp_spn_.365
## grp_sph.000
## grp_sp_.000
## grp_spc_.007
## grp_spn_.007
## grp_spc_.090
## grp_spn_.090
## grp_spc_.365 0.420
## grp_spn_.365 0.634 0.420
## grp_spct_r. 0.390 0.432 0.3903.7.2 Pairwise comparison
3.7.2.1 List of contrasts
List of contrast names and left and right coordinates on the X axis
Open code
contrast <- c(
## healthy vs. others
"healthy vs. no rejection D_0",
"healthy vs. no rejection D_7",
"healthy vs. no rejection D_90",
"healthy vs. no rejection D_365",
"healthy vs. acute rejection D_0",
"healthy vs. acute rejection D_7",
"healthy vs. acute rejection D_90",
"healthy vs. acute rejection D_365",
"healthy vs. acute rejection D_biopsy",
## between-subjects comparisons
"no rejection D_0 vs. acute rejection D_0",
"no rejection D_7 vs. acute rejection D_7",
"no rejection D_90 vs. acute rejection D_90",
"no rejection D_365 vs. acute rejection D_365",
"no rejection D_0 vs. acute rejection D_biopsy",
"no rejection D_7 vs. acute rejection D_biopsy",
"no rejection D_90 vs. acute rejection D_biopsy",
"no rejection D_365 vs. acute rejection D_biopsy",
## within-subjects comparisons
### no rejection
"no rejection D_0 vs. no rejection D_7",
"no rejection D_0 vs. no rejection D_90",
"no rejection D_0 vs. no rejection D_365",
"no rejection D_7 vs. no rejection D_90",
"no rejection D_7 vs. no rejection D_365",
"no rejection D_90 vs. no rejection D_365",
### acute rejection
"acute rejection D_0 vs. acute rejection D_7",
"acute rejection D_0 vs. acute rejection D_90",
"acute rejection D_0 vs. acute rejection D_365",
"acute rejection D_7 vs. acute rejection D_90",
"acute rejection D_7 vs. acute rejection D_365",
"acute rejection D_90 vs. acute rejection D_365",
"acute rejection D_0 vs. acute rejection D_biopsy",
"acute rejection D_7 vs. acute rejection D_biopsy",
"acute rejection D_90 vs. acute rejection D_biopsy",
"acute rejection D_365 vs. acute rejection D_biopsy"
)
xcoord_1 <- c(rep(6, 9), seq(0.8, 3.8, 1), seq(0.8, 3.8, 1),
rep(0.8, 3), rep(1.8, 2), 2.8,
rep(1.2, 3), rep(2.2, 2), 3.2,
seq(1.2, 4.2, 1))
xcoord_2 <- c(seq(0.8, 3.8, 1), seq(1.2, 4.2, 1),
5, seq(1.2, 4.2, 1), rep(5, 4),
seq(1.8, 3.8, 1), 2.8, 3.8, 3.8,
2.2, 3.2, 4.2, 3.2, 4.2, 4.2,
rep(5, 4))
anot_coord <- (xcoord_1+xcoord_2)/2
pair_table <- data.frame(contrast, xcoord_1, xcoord_2,
anot_coord)
pair_table <- pair_table %>%
mutate(term1 = if_else(xcoord_1 == 6, 'healthy.000', 'other')) %>%
mutate(term1 = if_else(xcoord_1 == 0.8, 'no_rejection.000', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 1.8, 'no_rejection.007', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 2.8, 'no_rejection.090', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 3.8, 'no_rejection.365', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 5, 'acute_rejection.biopsy', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 1.2, 'acute_rejection.000', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 2.2, 'acute_rejection.007', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 3.2, 'acute_rejection.090', term1)) %>%
mutate(term1 = if_else(xcoord_1 == 4.2, 'acute_rejection.365', term1)) %>%
mutate(term2 = if_else(xcoord_2 == 6, 'healthy.000', 'other')) %>%
mutate(term2 = if_else(xcoord_2 == 0.8, 'no_rejection.000', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 1.8, 'no_rejection.007', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 2.8, 'no_rejection.090', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 3.8, 'no_rejection.365', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 5, 'acute_rejection.biopsy', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 1.2, 'acute_rejection.000', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 2.2, 'acute_rejection.007', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 3.2, 'acute_rejection.090', term2)) %>%
mutate(term2 = if_else(xcoord_2 == 4.2, 'acute_rejection.365', term2))
pair_table$p <- NA3.7.2.2 Between-subjects comparisons
Open code
set.seed(16)
for (i in 1:17){
dat_pair <- model_data %>%
filter(group_sep == paste(pair_table$term1[i])|
group_sep == paste(pair_table$term2[i]))
g1 <- dat_pair %>% filter(group_sep == paste(pair_table$term1[i]))
g2 <- dat_pair %>% filter(group_sep == paste(pair_table$term2[i]))
pair_table$p[i] <- wilcox.test(g1$il36_b_value, g2$il36_b_value, paired=FALSE)$p.value
}3.7.2.3 Within-subject comparisons
Open code
set.seed(16)
for (i in 18:33){
dat_pair <- model_data %>%
filter(group_sep == paste(pair_table$term1[i]) |
group_sep == paste(pair_table$term2[i])) %>%
select(id, group_sep, il36_b_value) %>%
pivot_wider(names_from = group_sep, values_from = il36_b_value) %>%
filter(complete.cases(.)) %>% data.frame()
g1 <- dat_pair[,2]
g2 <- dat_pair[,3]
pair_table$p[i] <- wilcox.test(c(g1), c(g2), paired=TRUE)$p.value
}
pair_table <- pair_table[-c(14:15, 19:20, 22, 25:26, 28, 30:33),]
pair_table$cole <- c(rep("black", 15), rep("blue3", 3), rep("red4", 3))
pair_table$fdr <- p.adjust(pair_table$p, method = 'BH')3.7.2.4 Resulting table
Open code
## the whole `pair table`
pair_table
## contrast xcoord_1 xcoord_2 anot_coord
## 1 healthy vs. no rejection D_0 6.0 0.8 3.4
## 2 healthy vs. no rejection D_7 6.0 1.8 3.9
## 3 healthy vs. no rejection D_90 6.0 2.8 4.4
## 4 healthy vs. no rejection D_365 6.0 3.8 4.9
## 5 healthy vs. acute rejection D_0 6.0 1.2 3.6
## 6 healthy vs. acute rejection D_7 6.0 2.2 4.1
## 7 healthy vs. acute rejection D_90 6.0 3.2 4.6
## 8 healthy vs. acute rejection D_365 6.0 4.2 5.1
## 9 healthy vs. acute rejection D_biopsy 6.0 5.0 5.5
## 10 no rejection D_0 vs. acute rejection D_0 0.8 1.2 1.0
## 11 no rejection D_7 vs. acute rejection D_7 1.8 2.2 2.0
## 12 no rejection D_90 vs. acute rejection D_90 2.8 3.2 3.0
## 13 no rejection D_365 vs. acute rejection D_365 3.8 4.2 4.0
## 16 no rejection D_90 vs. acute rejection D_biopsy 2.8 5.0 3.9
## 17 no rejection D_365 vs. acute rejection D_biopsy 3.8 5.0 4.4
## 18 no rejection D_0 vs. no rejection D_7 0.8 1.8 1.3
## 21 no rejection D_7 vs. no rejection D_90 1.8 2.8 2.3
## 23 no rejection D_90 vs. no rejection D_365 2.8 3.8 3.3
## 24 acute rejection D_0 vs. acute rejection D_7 1.2 2.2 1.7
## 27 acute rejection D_7 vs. acute rejection D_90 2.2 3.2 2.7
## 29 acute rejection D_90 vs. acute rejection D_365 3.2 4.2 3.7
## term1 term2 p cole fdr
## 1 healthy.000 no_rejection.000 1.507542e-15 black 3.165839e-14
## 2 healthy.000 no_rejection.007 1.250083e-13 black 1.312587e-12
## 3 healthy.000 no_rejection.090 5.134600e-09 black 2.695665e-08
## 4 healthy.000 no_rejection.365 7.018490e-08 black 2.947766e-07
## 5 healthy.000 acute_rejection.000 1.419560e-11 black 9.936919e-11
## 6 healthy.000 acute_rejection.007 1.164937e-06 black 4.077280e-06
## 7 healthy.000 acute_rejection.090 3.199889e-02 black 4.799833e-02
## 8 healthy.000 acute_rejection.365 6.882559e-01 black 6.882559e-01
## 9 healthy.000 acute_rejection.biopsy 3.350608e-01 black 3.518138e-01
## 10 no_rejection.000 acute_rejection.000 1.489642e-01 black 1.747311e-01
## 11 no_rejection.007 acute_rejection.007 2.290071e-02 black 4.007624e-02
## 12 no_rejection.090 acute_rejection.090 8.967787e-04 black 2.354044e-03
## 13 no_rejection.365 acute_rejection.365 1.013818e-03 black 2.365575e-03
## 16 no_rejection.090 acute_rejection.biopsy 5.687729e-03 black 1.085839e-02
## 17 no_rejection.365 acute_rejection.biopsy 2.736817e-02 black 4.421013e-02
## 18 no_rejection.000 no_rejection.007 1.152029e-05 blue3 3.456087e-05
## 21 no_rejection.007 no_rejection.090 3.593050e-02 blue3 5.030269e-02
## 23 no_rejection.090 no_rejection.365 1.755844e-01 blue3 1.940670e-01
## 24 acute_rejection.000 acute_rejection.007 1.390606e-03 red4 2.920273e-03
## 27 acute_rejection.007 acute_rejection.090 5.156805e-02 red4 6.768306e-02
## 29 acute_rejection.090 acute_rejection.365 1.497695e-01 red4 1.747311e-01
## only significant differences
pair_table %>%
filter(p<0.05) %>%
select(contrast, p, fdr) %>%
mutate(star = if_else(p<0.05, '*', '')) %>%
mutate(star = if_else(p<0.01, '**', star)) %>%
mutate(star = if_else(p<0.001, '***', star)) %>%
data.frame()
## contrast p fdr
## 1 healthy vs. no rejection D_0 1.507542e-15 3.165839e-14
## 2 healthy vs. no rejection D_7 1.250083e-13 1.312587e-12
## 3 healthy vs. no rejection D_90 5.134600e-09 2.695665e-08
## 4 healthy vs. no rejection D_365 7.018490e-08 2.947766e-07
## 5 healthy vs. acute rejection D_0 1.419560e-11 9.936919e-11
## 6 healthy vs. acute rejection D_7 1.164937e-06 4.077280e-06
## 7 healthy vs. acute rejection D_90 3.199889e-02 4.799833e-02
## 8 no rejection D_7 vs. acute rejection D_7 2.290071e-02 4.007624e-02
## 9 no rejection D_90 vs. acute rejection D_90 8.967787e-04 2.354044e-03
## 10 no rejection D_365 vs. acute rejection D_365 1.013818e-03 2.365575e-03
## 11 no rejection D_90 vs. acute rejection D_biopsy 5.687729e-03 1.085839e-02
## 12 no rejection D_365 vs. acute rejection D_biopsy 2.736817e-02 4.421013e-02
## 13 no rejection D_0 vs. no rejection D_7 1.152029e-05 3.456087e-05
## 14 no rejection D_7 vs. no rejection D_90 3.593050e-02 5.030269e-02
## 15 acute rejection D_0 vs. acute rejection D_7 1.390606e-03 2.920273e-03
## star
## 1 ***
## 2 ***
## 3 ***
## 4 ***
## 5 ***
## 6 ***
## 7 *
## 8 *
## 9 ***
## 10 **
## 11 **
## 12 *
## 13 ***
## 14 *
## 15 **
pair_rest <- pair_table %>%
filter(p<0.05,
term1 != 'healthy.000') %>%
mutate(star = if_else(p<0.05, '*', '')) %>%
mutate(star = if_else(p<0.01, '**', star)) %>%
mutate(star = if_else(p<0.001, '***', star)) %>%
data.frame()Healthy does differ from all but acute rejection during biopsy and D365
3.7.3 Visualization
Open code
## plot
plotos <- data_long %>%
filter(!is.na(il36_b_value)) %>%
mutate(group = factor(group, levels =
c("no_rejection",
"acute_rejection",
"healthy")),
time_point = factor(time_point)) %>%
mutate(time_point = if_else(group == 'healthy', 'healthy', time_point)) %>%
mutate(group_sep = interaction(group, time_point)) %>%
## plotting
ggplot(aes(x=time_point, y=il36_b_value)) +
geom_boxplot(aes(fill=group, group = group_sep), outlier.shape = NA,
size=0.5,position = position_dodge(width = 0.84), color='black') +
labs(x = "Time post TX", y = expression("IL-36" * beta * " (pg/ml)")) +
scale_x_discrete(labels = c("Before TX", "7 days", "3 months", "365 days",
"Biopsy", "Healthy")) +
theme(axis.text = element_text(size = 11),
axis.title = element_text(size = 12),
legend.text = element_text(size = 11),
legend.title = element_text(size = 12)) +
scale_fill_manual(values = c("no_rejection" = "skyblue2",
"acute_rejection" = "coral",
"healthy" = "darkolivegreen3"),
labels = c("no rejection", "acute rejection", "healthy")) +
coord_cartesian(ylim = c(0, 15.5)) +
annotate("text", x = 6, y = 5, label = "#", size = 7)
pair_rest$y <- c(9.9, 7, 6.4, 9,
8.1, 15.5, 12.5, 11.5)
plotos <- plotos + geom_segment(data = pair_rest,
aes(x = xcoord_1, xend = xcoord_2, y = y, yend = y),
color = pair_rest$cole,
linetype = "solid", size = 0.4)
plot_il36_b <- plotos + annotate("text", x = pair_rest$anot_coord, y = pair_rest$y+0.15,
label = pair_rest$star, size = 7, color=pair_rest$cole)
plot_il36_b
4 Correlation of ILs and renal function markers
4.1 Simple correlation matrix based on patient data
Open code
tr_cor <- patient_table %>%
filter(group != 'healthy') %>%
select(egfr_mean:il36_b_mean)
cor.mtest(tr_cor, conf.level = 0.95)
## $p
## egfr_mean creatinine_mean il1_a_mean il1_b_mean
## egfr_mean 0.000000e+00 2.103729e-16 6.139177e-01 2.941143e-01
## creatinine_mean 2.103729e-16 0.000000e+00 9.054514e-01 9.825541e-01
## il1_a_mean 6.139177e-01 9.054514e-01 0.000000e+00 5.817931e-18
## il1_b_mean 2.941143e-01 9.825541e-01 5.817931e-18 0.000000e+00
## il1_ra_mean 7.275104e-01 8.430982e-01 4.288839e-15 2.316736e-11
## il18_mean 8.565717e-01 7.798943e-01 1.680311e-15 7.659848e-14
## il18_bp_mean 1.795693e-01 2.157285e-01 3.841527e-01 7.085025e-01
## il18_free_mean 6.704398e-01 6.534084e-01 3.949705e-15 1.055580e-11
## il36_b_mean 2.452979e-01 3.918567e-01 9.302473e-05 8.379317e-06
## il1_ra_mean il18_mean il18_bp_mean il18_free_mean
## egfr_mean 7.275104e-01 8.565717e-01 0.1795693 6.704398e-01
## creatinine_mean 8.430982e-01 7.798943e-01 0.2157285 6.534084e-01
## il1_a_mean 4.288839e-15 1.680311e-15 0.3841527 3.949705e-15
## il1_b_mean 2.316736e-11 7.659848e-14 0.7085025 1.055580e-11
## il1_ra_mean 0.000000e+00 2.176954e-20 0.5645638 2.259573e-13
## il18_mean 2.176954e-20 0.000000e+00 0.9768424 6.833531e-85
## il18_bp_mean 5.645638e-01 9.768424e-01 0.0000000 4.524752e-01
## il18_free_mean 2.259573e-13 6.833531e-85 0.4524752 0.000000e+00
## il36_b_mean 2.540206e-02 1.512619e-03 0.3023744 1.849524e-03
## il36_b_mean
## egfr_mean 2.452979e-01
## creatinine_mean 3.918567e-01
## il1_a_mean 9.302473e-05
## il1_b_mean 8.379317e-06
## il1_ra_mean 2.540206e-02
## il18_mean 1.512619e-03
## il18_bp_mean 3.023744e-01
## il18_free_mean 1.849524e-03
## il36_b_mean 0.000000e+00
##
## $lowCI
## egfr_mean creatinine_mean il1_a_mean il1_b_mean il1_ra_mean
## egfr_mean 1.00000000 -0.71827559 -0.1246886 -0.07833667 -0.19605109
## creatinine_mean -0.71827559 1.00000000 -0.1770066 -0.16527820 -0.18358576
## il1_a_mean -0.12468860 -0.17700661 1.0000000 0.54243007 0.48621354
## il1_b_mean -0.07833667 -0.16527820 0.5424301 1.00000000 0.39813645
## il1_ra_mean -0.19605109 -0.18358576 0.4862135 0.39813645 1.00000000
## il18_mean -0.18215813 -0.14367730 0.4947200 0.45880699 0.58405262
## il18_bp_mean -0.27670352 -0.06215646 -0.2387838 -0.19815124 -0.21478835
## il18_free_mean -0.20241378 -0.12937782 0.4869689 0.40705248 0.44800513
## il36_b_mean -0.06868499 -0.23764656 0.1685782 0.21530828 0.02392036
## il18_mean il18_bp_mean il18_free_mean il36_b_mean
## egfr_mean -0.1821581 -0.27670352 -0.2024138 -0.06868499
## creatinine_mean -0.1436773 -0.06215646 -0.1293778 -0.23764656
## il1_a_mean 0.4947200 -0.23878380 0.4869689 0.16857815
## il1_b_mean 0.4588070 -0.19815124 0.4070525 0.21530828
## il1_ra_mean 0.5840526 -0.21478835 0.4480051 0.02392036
## il18_mean 1.0000000 -0.16467977 0.9574656 0.10512859
## il18_bp_mean -0.1646798 1.00000000 -0.2291107 -0.07985789
## il18_free_mean 0.9574656 -0.22911074 1.0000000 0.10003378
## il36_b_mean 0.1051286 -0.07985789 0.1000338 1.00000000
##
## $uppCI
## egfr_mean creatinine_mean il1_a_mean il1_b_mean il1_ra_mean
## egfr_mean 1.00000000 -0.51293155 0.20891140 0.2532441 0.1378676
## creatinine_mean -0.51293155 1.00000000 0.15716924 0.1689303 0.1505300
## il1_a_mean 0.20891140 0.15716924 1.00000000 0.7374956 0.7005925
## il1_b_mean 0.25324409 0.16893033 0.73749557 1.0000000 0.6403732
## il1_ra_mean 0.13786762 0.15052999 0.70059252 0.6403732 1.0000000
## il18_mean 0.15197323 0.19034518 0.70625117 0.6821760 0.7640880
## il18_bp_mean 0.05321464 0.26839831 0.09361479 0.1357235 0.1186266
## il18_free_mean 0.13136199 0.20434877 0.70109612 0.6466077 0.6748387
## il36_b_mean 0.26230338 0.09480974 0.46805481 0.5050816 0.3463574
## il18_mean il18_bp_mean il18_free_mean il36_b_mean
## egfr_mean 0.1519732 0.05321464 0.1313620 0.26230338
## creatinine_mean 0.1903452 0.26839831 0.2043488 0.09480974
## il1_a_mean 0.7062512 0.09361479 0.7010961 0.46805481
## il1_b_mean 0.6821760 0.13572347 0.6466077 0.50508160
## il1_ra_mean 0.7640880 0.11862657 0.6748387 0.34635737
## il18_mean 1.0000000 0.16952788 0.9781094 0.41603888
## il18_bp_mean 0.1695279 1.00000000 0.1037483 0.25181091
## il18_free_mean 0.9781094 0.10374828 1.0000000 0.41177200
## il36_b_mean 0.4160389 0.25181091 0.4117720 1.00000000
tr_cor <- cor(tr_cor, method='spearman',
use = 'pairwise.complete.obs')
corrplot(tr_cor, addCoef.col = 'black',tl.col = "white")
text(1:9, 9.9, expression("eGFR", "Creatinine", "IL-1" * alpha,
"IL-1" * beta, "IL-1 RA", "IL-18", "IL-18 BP",
"IL-18 (free)", "IL-36" * beta))
text(-0.1, 9:1, expression("eGFR", "Creatinine", "IL-1" * alpha,
"IL-1" * beta, "IL-1 RA", "IL-18", "IL-18 BP",
"IL-18 (free)", "IL-1" * beta))
There is not any significant correlation between interleukins and markers of renal function. At least modest (but still not statistically significant) correlation is for interleukin-18BP.
4.2 Adjsuting IL-18BP model for renal function
To further validate that renal functions does not play an important role for above-shown findings, we will re-fit mixed-effects model adjusting also for the effect of renal function. We will do it only for IL-18BP where modest correlation between average values of renal markers and IL-18BP levels was found.
4.2.1 Adjustment for EGFR
Open code
model_data <- data_long %>%
filter(!is.na(il18_bp_value))
## Complex models (incuding healthy and biopsy data)
model_il18_bp <- lmer(log10(il18_bp_value) ~
group_sep +
(1|id), data = model_data, REML = FALSE)
model_il18_bp_adjEGFR <- lmer(log10(il18_bp_value) ~ log2(egfr) +
group_sep +
(1|id), data = model_data, REML = FALSE)
model_il18_bp_adjEGFR_REML <- lmer(log10(il18_bp_value) ~ log2(egfr) +
group_sep +
(1|id), data = model_data, REML = TRUE)
Anova(model_il18_bp_adjEGFR_REML, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il18_bp_value)
## F Df Df.res Pr(>F)
## log2(egfr) 7.9151 1 547.11 0.005079 **
## group_sep 0.8430 9 498.72 0.576683
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_il18_bp_adjEGFR_REML)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il18_bp_value) ~ log2(egfr) + group_sep + (1 | id)
## Data: model_data
##
## REML criterion at convergence: -484.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -10.0627 -0.3090 0.0389 0.3262 4.7340
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.00581 0.07622
## Residual 0.01870 0.13675
## Number of obs: 577, groups: id, 169
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.541856 0.035691 99.236
## log2(egfr) -0.023511 0.008321 -2.825
## group_sephealthy.000 -0.082453 0.047640 -1.731
## group_sepno_rejection.000 0.003293 0.029488 0.112
## group_sepacute_rejection.007 -0.022998 0.032948 -0.698
## group_sepno_rejection.007 -0.008345 0.032425 -0.257
## group_sepacute_rejection.090 -0.061574 0.037176 -1.656
## group_sepno_rejection.090 -0.043781 0.037677 -1.162
## group_sepacute_rejection.365 -0.030476 0.038551 -0.791
## group_sepno_rejection.365 -0.054496 0.038220 -1.426
## group_sepacute_rejection.biopsy -0.015365 0.036709 -0.419
##
## Correlation of Fixed Effects:
## (Intr) lg2(g) g_.000 g__.000 grp_spc_.007 grp_spn_.007
## log2(egfr) 0.724
## grp_sph.000 -0.805 -0.620
## grp_sp_.000 -0.568 0.011 0.425
## grp_spc_.007 -0.703 -0.414 0.559 0.484
## grp_spn_.007 -0.827 -0.419 0.652 0.694 0.617
## grp_spc_.090 -0.769 -0.568 0.620 0.427 0.623 0.632
## grp_spn_.090 -0.895 -0.613 0.718 0.595 0.636 0.804
## grp_spc_.365 -0.755 -0.566 0.610 0.412 0.609 0.618
## grp_spn_.365 -0.894 -0.621 0.718 0.586 0.634 0.799
## grp_spct_r. -0.622 -0.364 0.494 0.430 0.539 0.547
## grp_spc_.090 grp_spn_.090 grp_spc_.365 grp_spn_.365
## log2(egfr)
## grp_sph.000
## grp_sp_.000
## grp_spc_.007
## grp_spn_.007
## grp_spc_.090
## grp_spn_.090 0.687
## grp_spc_.365 0.653 0.675
## grp_spn_.365 0.687 0.845 0.674
## grp_spct_r. 0.551 0.563 0.535 0.561
AIC(model_il18_bp, model_il18_bp_adjEGFR)
## df AIC
## model_il18_bp 12 -519.0147
## model_il18_bp_adjEGFR 13 -525.0759
BIC(model_il18_bp, model_il18_bp_adjEGFR)
## df BIC
## model_il18_bp 12 -466.7206
## model_il18_bp_adjEGFR 13 -468.4240
## Patient-only, planned-measures-only models
model_data <- data_long %>%
filter(!is.na(il18_bp_value),
group != 'healthy',
time_point != 'biopsy')
model_il18_bp <- lmer(log10(il18_bp_value) ~
group*time_point +
(1|id), data = model_data, REML = FALSE)
model_il18_bp_adjEGFR <- lmer(log10(il18_bp_value) ~ log2(egfr) +
group*time_point +
(1|id), data = model_data, REML = FALSE)
model_il18_bp_adjEGFR_REML <- lmer(log10(il18_bp_value) ~ log2(egfr) +
group*time_point +
(1|id), data = model_data, REML = TRUE)
Anova(model_il18_bp_adjEGFR_REML, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il18_bp_value)
## F Df Df.res Pr(>F)
## log2(egfr) 8.3333 1 493.34 0.004063 **
## group 0.0399 1 136.84 0.841907
## time_point 1.1440 3 427.71 0.330982
## group:time_point 0.4743 3 383.95 0.700327
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_il18_bp_adjEGFR_REML)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il18_bp_value) ~ log2(egfr) + group * time_point + (1 | id)
## Data: model_data
##
## REML criterion at convergence: -395.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.6489 -0.3011 0.0298 0.3526 4.5301
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.006333 0.07958
## Residual 0.020298 0.14247
## Number of obs: 518, groups: id, 138
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.537355 0.032864 107.637
## log2(egfr) -0.025990 0.008963 -2.900
## groupacute_rejection -0.003096 0.030741 -0.101
## time_point007 -0.007498 0.025666 -0.292
## time_point090 -0.040094 0.033041 -1.213
## time_point365 -0.050629 0.033752 -1.500
## groupacute_rejection:time_point007 -0.011536 0.037562 -0.307
## groupacute_rejection:time_point090 -0.015374 0.038446 -0.400
## groupacute_rejection:time_point365 0.027516 0.039657 0.694
##
## Correlation of Fixed Effects:
## (Intr) lg2(g) grpct_ tm_007 tm_090 tm_365 g_:_00 g_:_09
## log2(egfr) 0.857
## grpct_rjctn -0.293 -0.012
## time_pnt007 -0.759 -0.583 0.284
## time_pnt090 -0.856 -0.764 0.224 0.703
## time_pnt365 -0.855 -0.768 0.220 0.700 0.783
## grpct_:_007 0.184 0.008 -0.635 -0.455 -0.182 -0.178
## grpct_:_090 0.230 0.066 -0.622 -0.260 -0.409 -0.220 0.509
## grpct_:_365 0.218 0.058 -0.604 -0.249 -0.212 -0.395 0.494 0.486
AIC(model_il18_bp, model_il18_bp_adjEGFR)
## df AIC
## model_il18_bp 10 -422.7043
## model_il18_bp_adjEGFR 11 -429.1605
BIC(model_il18_bp, model_il18_bp_adjEGFR)
## df BIC
## model_il18_bp 10 -380.2045
## model_il18_bp_adjEGFR 11 -382.4108Inclusion of EGFR improves the fit noticeably
4.2.2 Adjustment for creatinine
Open code
## Complex models (including healthy and biopsy data)
model_data <- data_long %>%
filter(!is.na(il18_bp_value))
model_il18_bp <- lmer(log10(il18_bp_value) ~
group_sep +
(1|id), data = model_data, REML = FALSE)
model_il18_bp_adjCRE <- lmer(log2(il18_bp_value) ~ log2(creatinine) +
group_sep +
(1|id), data = model_data, REML = FALSE)
model_il18_bp_adjCRE_REML <- lmer(log10(il18_bp_value) ~ log2(creatinine) +
group_sep +
(1|id), data = model_data, REML = TRUE)
Anova(model_il18_bp_adjCRE_REML, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il18_bp_value)
## F Df Df.res Pr(>F)
## log2(creatinine) 8.4351 1 534.24 0.003833 **
## group_sep 0.8406 9 499.16 0.578871
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_il18_bp_adjCRE_REML)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il18_bp_value) ~ log2(creatinine) + group_sep + (1 | id)
## Data: model_data
##
## REML criterion at convergence: -485.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -10.0607 -0.3099 0.0346 0.3382 4.7544
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.005777 0.0760
## Residual 0.018700 0.1367
## Number of obs: 577, groups: id, 169
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.345720 0.095493 35.036
## log2(creatinine) 0.028785 0.009868 2.917
## group_sephealthy.000 -0.080801 0.047398 -1.705
## group_sepno_rejection.000 0.004245 0.029466 0.144
## group_sepacute_rejection.007 -0.022510 0.032840 -0.685
## group_sepno_rejection.007 -0.006920 0.032434 -0.213
## group_sepacute_rejection.090 -0.060569 0.036995 -1.637
## group_sepno_rejection.090 -0.041914 0.037610 -1.114
## group_sepacute_rejection.365 -0.029155 0.038422 -0.759
## group_sepno_rejection.365 -0.052324 0.038205 -1.370
## group_sepacute_rejection.biopsy -0.014570 0.036653 -0.398
##
## Correlation of Fixed Effects:
## (Intr) lg2(c) g_.000 g__.000 grp_spc_.007 grp_spn_.007
## log2(crtnn) -0.966
## grp_sph.000 -0.729 0.616
## grp_sp_.000 -0.216 0.001 0.434
## grp_spc_.007 -0.545 0.408 0.556 0.490
## grp_spn_.007 -0.602 0.421 0.653 0.698 0.617
## grp_spc_.090 -0.678 0.562 0.617 0.436 0.620 0.632
## grp_spn_.090 -0.760 0.612 0.716 0.602 0.633 0.804
## grp_spc_.365 -0.672 0.562 0.607 0.421 0.606 0.618
## grp_spn_.365 -0.766 0.621 0.717 0.593 0.631 0.799
## grp_spct_r. -0.483 0.360 0.492 0.435 0.537 0.547
## grp_spc_.090 grp_spn_.090 grp_spc_.365 grp_spn_.365
## log2(crtnn)
## grp_sph.000
## grp_sp_.000
## grp_spc_.007
## grp_spn_.007
## grp_spc_.090
## grp_spn_.090 0.685
## grp_spc_.365 0.650 0.673
## grp_spn_.365 0.685 0.845 0.673
## grp_spct_r. 0.549 0.561 0.533 0.559
AIC(model_il18_bp, model_il18_bp_adjCRE)
## df AIC
## model_il18_bp 12 -519.0147
## model_il18_bp_adjCRE 13 859.8325
BIC(model_il18_bp, model_il18_bp_adjCRE)
## df BIC
## model_il18_bp 12 -466.7206
## model_il18_bp_adjCRE 13 916.4844
## Patient-only, planned-measures-only models
model_data <- data_long %>%
filter(!is.na(il18_bp_value),
group != 'healthy',
time_point != 'biopsy')
model_il18_bp <- lmer(log10(il18_bp_value) ~
group*time_point +
(1|id), data = model_data, REML = FALSE)
model_il18_bp_adjEGFR <- lmer(log10(il18_bp_value) ~ log2(egfr) +
group*time_point +
(1|id), data = model_data, REML = FALSE)
model_il18_bp_adjEGFR_REML <- lmer(log10(il18_bp_value) ~ log2(egfr) +
group*time_point +
(1|id), data = model_data, REML = TRUE)
Anova(model_il18_bp_adjEGFR_REML, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il18_bp_value)
## F Df Df.res Pr(>F)
## log2(egfr) 8.3333 1 493.34 0.004063 **
## group 0.0399 1 136.84 0.841907
## time_point 1.1440 3 427.71 0.330982
## group:time_point 0.4743 3 383.95 0.700327
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(model_il18_bp, model_il18_bp_adjEGFR)
## df AIC
## model_il18_bp 10 -422.7043
## model_il18_bp_adjEGFR 11 -429.1605
BIC(model_il18_bp, model_il18_bp_adjEGFR)
## df BIC
## model_il18_bp 10 -380.2045
## model_il18_bp_adjEGFR 11 -382.4108The analysis shows that when renal function are accounted for (included as a covariate to the mixed-effects model), there is no remaining variability explainable by other factors such as group [acute rejection/no rejection/healthy], time_point and group:timepoint interaction. It suggests that IL-18BP might be driven mainly with renal function and its change over time whereas the effects of group and timepoint etc may be indirect.
4.3 Adjusted free IL-18 model
Given that IL-18BP is affected with renal function, and the IL18-BP affects level of other forms of IL-18 (including free IL-18), we will fit model explaining the IL-18 levels but simultaneously adjusting for renal functions
4.3.1 Adjustment for EGFR
Open code
model_data <- data_long %>%
filter(!is.na(il18_free))
## Complex models (incuding healthy and biopsy data)
model_il18_free <- lmer(log10(il18_free) ~
group_sep +
(1|id), data = model_data, REML = FALSE)
model_il18_free_adjEGFR <- lmer(log10(il18_free) ~ log2(egfr) +
group_sep +
(1|id), data = model_data, REML = FALSE)
model_il18_free_adjEGFR_REML <- lmer(log10(il18_free) ~ log2(egfr) +
group_sep +
(1|id), data = model_data, REML = TRUE)
Anova(model_il18_free_adjEGFR_REML, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il18_free)
## F Df Df.res Pr(>F)
## log2(egfr) 0.6484 1 454.0 0.4211
## group_sep 18.2407 9 408.3 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(model_il18_free, model_il18_free_adjEGFR)
## df AIC
## model_il18_free 12 227.9601
## model_il18_free_adjEGFR 13 229.2903
BIC(model_il18_free, model_il18_free_adjEGFR)
## df BIC
## model_il18_free 12 278.0706
## model_il18_free_adjEGFR 13 283.5766
## contrasts
em_il18_free <- emmeans(model_il18_free_adjEGFR_REML, specs = pairwise ~ group_sep,
adjust = 'none', type = 'response')
em_il18_free_cont <- data.frame(em_il18_free$contrasts)[c(
1:3, 10:17, 19, 25, 26, 32, 36, 37, 39, 42, 43, 45
),]
as_tibble(em_il18_free_cont) %>%
filter(p.value<0.05) %>% data.frame()
## contrast ratio SE df null
## 1 acute_rejection.000 / healthy.000 2.4073392 0.53767731 434.9662 1
## 2 acute_rejection.000 / no_rejection.000 0.6793931 0.09351587 425.4130 1
## 3 acute_rejection.000 / acute_rejection.007 1.7959185 0.27628214 368.5469 1
## 4 healthy.000 / no_rejection.000 0.2822174 0.05813244 432.7976 1
## 5 healthy.000 / no_rejection.007 0.6344966 0.10611335 418.2433 1
## 6 healthy.000 / acute_rejection.090 0.3025239 0.05310212 416.0604 1
## 7 healthy.000 / no_rejection.090 0.3225245 0.05607866 443.6156 1
## 8 healthy.000 / acute_rejection.365 0.3762615 0.06780092 424.0318 1
## 9 healthy.000 / no_rejection.365 0.4349602 0.07537710 443.7478 1
## 10 healthy.000 / acute_rejection.biopsy 0.3113010 0.06239168 447.3936 1
## 11 no_rejection.000 / no_rejection.007 2.2482543 0.25255124 406.6198 1
## 12 acute_rejection.007 / acute_rejection.090 0.4055182 0.05689067 337.8121 1
## 13 no_rejection.007 / no_rejection.090 0.5083156 0.06651835 391.4777 1
## t.ratio p.value
## 1 3.933401 9.745268e-05
## 2 -2.808326 5.209278e-03
## 3 3.806037 1.653100e-04
## 4 -6.141612 1.851593e-09
## 5 -2.720179 6.796665e-03
## 6 -6.811330 3.398621e-11
## 7 -6.508020 2.061651e-10
## 8 -5.424479 9.787394e-08
## 9 -4.803908 2.132018e-06
## 10 -5.822679 1.106736e-08
## 11 7.212129 2.710851e-12
## 12 -6.433682 4.262068e-10
## 13 -5.170801 3.726822e-07
## Patient-only, planned-measures-only models
model_data <- data_long %>%
filter(!is.na(il18_free),
group != 'healthy',
time_point != 'biopsy')
model_il18_free <- lmer(log10(il18_free) ~
group*time_point +
(1|id), data = model_data, REML = FALSE)
model_il18_free_adjEGFR <- lmer(log10(il18_free) ~ log2(egfr) +
group*time_point +
(1|id), data = model_data, REML = FALSE)
model_il18_free_adjEGFR_REML <- lmer(log10(il18_free) ~ log2(egfr) +
group*time_point +
(1|id), data = model_data, REML = TRUE)
Anova(model_il18_free_adjEGFR_REML, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il18_free)
## F Df Df.res Pr(>F)
## log2(egfr) 1.1034 1 401.46 0.29415
## group 1.7354 1 127.35 0.19008
## time_point 38.0225 3 333.84 < 2e-16 ***
## group:time_point 2.9089 3 300.10 0.03485 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(model_il18_free, model_il18_free_adjEGFR)
## df AIC
## model_il18_free 10 228.2280
## model_il18_free_adjEGFR 11 229.0902
BIC(model_il18_free, model_il18_free_adjEGFR)
## df BIC
## model_il18_free 10 268.6781
## model_il18_free_adjEGFR 11 273.5853Inclusion of EGFR improves the fit noticeably
4.3.2 Adjustment for creatinine
Open code
## Complex models (including healthy and biopsy data)
model_data <- data_long %>%
filter(!is.na(il18_free))
model_il18_free <- lmer(log10(il18_free) ~
group_sep +
(1|id), data = model_data, REML = FALSE)
model_il18_free_adjCRE <- lmer(log10(il18_free) ~ log2(creatinine) +
group_sep +
(1|id), data = model_data, REML = FALSE)
model_il18_free_adjCRE_REML <- lmer(log10(il18_free) ~ log2(creatinine) +
group_sep +
(1|id), data = model_data, REML = TRUE)
Anova(model_il18_free_adjCRE_REML, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il18_free)
## F Df Df.res Pr(>F)
## log2(creatinine) 0.0606 1 445.85 0.8057
## group_sep 18.0057 9 408.47 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(model_il18_free, model_il18_free_adjCRE)
## df AIC
## model_il18_free 12 227.9601
## model_il18_free_adjCRE 13 229.8973
BIC(model_il18_free, model_il18_free_adjCRE)
## df BIC
## model_il18_free 12 278.0706
## model_il18_free_adjCRE 13 284.1836
## contrasts
em_il18_free <- emmeans(model_il18_free_adjCRE_REML, specs = pairwise ~ group_sep,
adjust = 'none', type = 'response')
em_il18_free_cont <- data.frame(em_il18_free$contrasts)[c(
1:3, 10:17, 19, 25, 26, 32, 36, 37, 39, 42, 43, 45
),]
as_tibble(em_il18_free_cont) %>%
filter(p.value<0.05) %>% data.frame()
## contrast ratio SE df null
## 1 acute_rejection.000 / healthy.000 2.2241926 0.49446636 431.4742 1
## 2 acute_rejection.000 / no_rejection.000 0.6801853 0.09368936 425.2481 1
## 3 acute_rejection.000 / acute_rejection.007 1.7300734 0.26527877 368.8704 1
## 4 healthy.000 / no_rejection.000 0.3058122 0.06224851 427.1593 1
## 5 healthy.000 / no_rejection.007 0.6623566 0.11018206 415.2399 1
## 6 healthy.000 / acute_rejection.090 0.3093849 0.05432380 414.9301 1
## 7 healthy.000 / no_rejection.090 0.3284553 0.05704034 442.9723 1
## 8 healthy.000 / acute_rejection.365 0.3842716 0.06922734 423.0664 1
## 9 healthy.000 / no_rejection.365 0.4421845 0.07652278 443.2303 1
## 10 healthy.000 / acute_rejection.biopsy 0.3250711 0.06504205 445.2994 1
## 11 no_rejection.000 / no_rejection.007 2.1658930 0.24151755 408.1660 1
## 12 acute_rejection.007 / acute_rejection.090 0.3977471 0.05580360 337.5963 1
## 13 no_rejection.007 / no_rejection.090 0.4958889 0.06485914 391.3846 1
## t.ratio p.value
## 1 3.595808 3.607667e-04
## 2 -2.797933 5.376774e-03
## 3 3.574970 3.968385e-04
## 4 -5.820564 1.152945e-08
## 5 -2.476434 1.366773e-02
## 6 -6.681433 7.635869e-11
## 7 -6.411028 3.708165e-10
## 8 -5.308879 1.784928e-07
## 9 -4.715392 3.236763e-06
## 10 -5.616153 3.443165e-08
## 11 6.930648 1.642603e-11
## 12 -6.571234 1.892147e-10
## 13 -5.362670 1.405310e-07
## Patient-only, planned-measures-only models
model_data <- data_long %>%
filter(!is.na(il18_free),
group != 'healthy',
time_point != 'biopsy')
model_il18_free <- lmer(log10(il18_free) ~
group*time_point +
(1|id), data = model_data, REML = FALSE)
model_il18_free_adjEGFR <- lmer(log10(il18_free) ~ log2(egfr) +
group*time_point +
(1|id), data = model_data, REML = FALSE)
model_il18_free_adjEGFR_REML <- lmer(log10(il18_free) ~ log2(egfr) +
group*time_point +
(1|id), data = model_data, REML = TRUE)
Anova(model_il18_free_adjEGFR_REML, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il18_free)
## F Df Df.res Pr(>F)
## log2(egfr) 1.1034 1 401.46 0.29415
## group 1.7354 1 127.35 0.19008
## time_point 38.0225 3 333.84 < 2e-16 ***
## group:time_point 2.9089 3 300.10 0.03485 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(model_il18_free, model_il18_free_adjEGFR)
## df AIC
## model_il18_free 10 228.2280
## model_il18_free_adjEGFR 11 229.0902
BIC(model_il18_free, model_il18_free_adjEGFR)
## df BIC
## model_il18_free 10 268.6781
## model_il18_free_adjEGFR 11 273.5853Adjustment for renal function does not modify result of the model explaining free IL-18. Markers of renal function do not clearly affect the pattern of difference in free IL-18 levels across groups or time points
What about total IL-18?
4.4 Adjusted total IL-18 model
4.4.1 Adjustment for EGFR
Open code
model_data <- data_long %>%
filter(!is.na(il18_value))
## Complex models (incuding healthy and biopsy data)
model_il18_value <- lmer(log10(il18_value) ~
group_sep +
(1|id), data = model_data, REML = FALSE)
model_il18_value_adjEGFR <- lmer(log10(il18_value) ~ log2(egfr) +
group_sep +
(1|id), data = model_data, REML = FALSE)
model_il18_value_adjEGFR_REML <- lmer(log10(il18_value) ~ log2(egfr) +
group_sep +
(1|id), data = model_data, REML = TRUE)
Anova(model_il18_value_adjEGFR_REML, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il18_value)
## F Df Df.res Pr(>F)
## log2(egfr) 0.0518 1 479.60 0.8201
## group_sep 20.2839 9 420.04 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(model_il18_value, model_il18_value_adjEGFR)
## df AIC
## model_il18_value 12 271.2115
## model_il18_value_adjEGFR 13 273.1579
BIC(model_il18_value, model_il18_value_adjEGFR)
## df BIC
## model_il18_value 12 321.7868
## model_il18_value_adjEGFR 13 327.9478
## contrasts
em_il18_value <- emmeans(model_il18_value_adjEGFR_REML, specs = pairwise ~ group_sep,
adjust = 'none', type = 'response')
em_il18_value_cont <- data.frame(em_il18_value$contrasts)[c(
1:3, 10:17, 19, 25, 26, 32, 36, 37, 39, 42, 43, 45
),]
as_tibble(em_il18_value_cont) %>%
filter(p.value<0.05) %>% data.frame()
## contrast ratio SE df null
## 1 acute_rejection.000 / healthy.000 2.8016661 0.60639470 454.4073 1
## 2 acute_rejection.000 / no_rejection.000 0.6587990 0.09349772 431.3921 1
## 3 acute_rejection.000 / acute_rejection.007 1.8091393 0.28340410 375.4448 1
## 4 healthy.000 / no_rejection.000 0.2351455 0.04639809 455.5534 1
## 5 healthy.000 / acute_rejection.007 0.6457369 0.11290118 430.7954 1
## 6 healthy.000 / no_rejection.007 0.5400838 0.08320735 434.4828 1
## 7 healthy.000 / acute_rejection.090 0.2686464 0.04407996 427.5350 1
## 8 healthy.000 / no_rejection.090 0.2750560 0.04449238 467.0760 1
## 9 healthy.000 / acute_rejection.365 0.3259036 0.05519896 439.1190 1
## 10 healthy.000 / no_rejection.365 0.3689904 0.05944705 467.1907 1
## 11 healthy.000 / acute_rejection.biopsy 0.2703141 0.05182941 468.8662 1
## 12 no_rejection.000 / no_rejection.007 2.2968076 0.26377620 414.3946 1
## 13 acute_rejection.007 / acute_rejection.090 0.4160308 0.05972821 342.3573 1
## 14 no_rejection.007 / no_rejection.090 0.5092838 0.06839126 396.4367 1
## t.ratio p.value
## 1 4.759798 2.608696e-06
## 2 -2.940618 3.451797e-03
## 3 3.784527 1.791956e-04
## 4 -7.336186 1.014798e-12
## 5 -2.501493 1.273714e-02
## 6 -3.998545 7.488620e-05
## 7 -8.010396 1.096085e-14
## 8 -7.979724 1.144495e-14
## 9 -6.619473 1.053554e-10
## 10 -6.188327 1.330858e-09
## 11 -6.822709 2.768522e-11
## 12 7.240387 2.192254e-12
## 13 -6.108626 2.730562e-09
## 14 -5.024607 7.647530e-07
## Patient-only, planned-measures-only models
model_data <- data_long %>%
filter(!is.na(il18_value),
group != 'healthy',
time_point != 'biopsy')
model_il18_value <- lmer(log10(il18_value) ~
group*time_point +
(1|id), data = model_data, REML = FALSE)
model_il18_value_adjEGFR <- lmer(log10(il18_value) ~ log2(egfr) +
group*time_point +
(1|id), data = model_data, REML = FALSE)
model_il18_value_adjEGFR_REML <- lmer(log10(il18_value) ~ log2(egfr) +
group*time_point +
(1|id), data = model_data, REML = TRUE)
Anova(model_il18_value_adjEGFR_REML, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il18_value)
## F Df Df.res Pr(>F)
## log2(egfr) 0.3786 1 404.44 0.53871
## group 2.4964 1 127.58 0.11658
## time_point 38.6155 3 335.53 < 2e-16 ***
## group:time_point 3.0134 3 301.04 0.03035 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(model_il18_value, model_il18_value_adjEGFR)
## df AIC
## model_il18_value 10 229.1316
## model_il18_value_adjEGFR 11 230.7408
BIC(model_il18_value, model_il18_value_adjEGFR)
## df BIC
## model_il18_value 10 269.6289
## model_il18_value_adjEGFR 11 275.28794.4.2 Adjustment for creatinine
Open code
## Complex models (including healthy and biopsy data)
model_data <- data_long %>%
filter(!is.na(il18_value))
model_il18_value <- lmer(log10(il18_value) ~
group_sep +
(1|id), data = model_data, REML = FALSE)
model_il18_value_adjCRE <- lmer(log10(il18_value) ~ log2(creatinine) +
group_sep +
(1|id), data = model_data, REML = FALSE)
model_il18_value_adjCRE_REML <- lmer(log10(il18_value) ~ log2(creatinine) +
group_sep +
(1|id), data = model_data, REML = TRUE)
Anova(model_il18_value_adjCRE_REML, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il18_value)
## F Df Df.res Pr(>F)
## log2(creatinine) 0.0663 1 471.87 0.7969
## group_sep 20.0253 9 420.33 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(model_il18_value, model_il18_value_adjCRE)
## df AIC
## model_il18_value 12 271.2115
## model_il18_value_adjCRE 13 273.1441
BIC(model_il18_value, model_il18_value_adjCRE)
## df BIC
## model_il18_value 12 321.7868
## model_il18_value_adjCRE 13 327.9341
## contrasts
em_il18_value <- emmeans(model_il18_value_adjCRE_REML, specs = pairwise ~ group_sep,
adjust = 'none', type = 'response')
em_il18_value_cont <- data.frame(em_il18_value$contrasts)[c(
1:3, 10:17, 19, 25, 26, 32, 36, 37, 39, 42, 43, 45
),]
as_tibble(em_il18_value_cont) %>%
filter(p.value<0.05) %>% data.frame()
## contrast ratio SE df null
## 1 acute_rejection.000 / healthy.000 2.6120988 0.56418523 450.2078 1
## 2 acute_rejection.000 / no_rejection.000 0.6587692 0.09347288 431.7167 1
## 3 acute_rejection.000 / acute_rejection.007 1.7506807 0.27327263 376.1026 1
## 4 healthy.000 / no_rejection.000 0.2521992 0.04927002 448.6505 1
## 5 healthy.000 / acute_rejection.007 0.6702199 0.11738276 429.1309 1
## 6 healthy.000 / no_rejection.007 0.5606891 0.08608655 430.9211 1
## 7 healthy.000 / acute_rejection.090 0.2740685 0.04503086 426.8454 1
## 8 healthy.000 / no_rejection.090 0.2794271 0.04516298 466.4659 1
## 9 healthy.000 / acute_rejection.365 0.3320203 0.05626098 438.3287 1
## 10 healthy.000 / no_rejection.365 0.3741927 0.06021292 466.7281 1
## 11 healthy.000 / acute_rejection.biopsy 0.2808052 0.05384462 466.7466 1
## 12 no_rejection.000 / no_rejection.007 2.2231996 0.25325311 416.6459 1
## 13 acute_rejection.007 / acute_rejection.090 0.4089233 0.05870204 342.3214 1
## 14 no_rejection.007 / no_rejection.090 0.4983638 0.06686930 396.7331 1
## t.ratio p.value
## 1 4.445379 1.106241e-05
## 2 -2.941585 3.441110e-03
## 3 3.587587 3.778455e-04
## 4 -7.051215 6.743530e-12
## 5 -2.284731 2.281589e-02
## 6 -3.768398 1.871701e-04
## 7 -7.877888 2.790387e-14
## 8 -7.888616 2.190487e-14
## 9 -6.506677 2.102148e-10
## 10 -6.108748 2.118180e-09
## 11 -6.623671 9.669253e-11
## 12 7.013614 9.438174e-12
## 13 -6.229265 1.374605e-09
## 14 -5.190319 3.358860e-07
## Patient-only, planned-measures-only models
model_data <- data_long %>%
filter(!is.na(il18_value),
group != 'healthy',
time_point != 'biopsy')
model_il18_value <- lmer(log10(il18_value) ~
group*time_point +
(1|id), data = model_data, REML = FALSE)
model_il18_value_adjEGFR <- lmer(log10(il18_value) ~ log2(egfr) +
group*time_point +
(1|id), data = model_data, REML = FALSE)
model_il18_value_adjEGFR_REML <- lmer(log10(il18_value) ~ log2(egfr) +
group*time_point +
(1|id), data = model_data, REML = TRUE)
Anova(model_il18_value_adjEGFR_REML, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il18_value)
## F Df Df.res Pr(>F)
## log2(egfr) 0.3786 1 404.44 0.53871
## group 2.4964 1 127.58 0.11658
## time_point 38.6155 3 335.53 < 2e-16 ***
## group:time_point 3.0134 3 301.04 0.03035 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(model_il18_value, model_il18_value_adjEGFR)
## df AIC
## model_il18_value 10 229.1316
## model_il18_value_adjEGFR 11 230.7408
BIC(model_il18_value, model_il18_value_adjEGFR)
## df BIC
## model_il18_value 10 269.6289
## model_il18_value_adjEGFR 11 275.2879Adjustment for renal function does not modify the results of the models explaining total IL-18. Markers of renal function do not clearly affect the above-shown patterns in total IL-18 levels and their difference across groups or time points
4.5 Renal function adjusment for other interleukins
4.5.1 IL-1a
Open code
model_data <- data_long %>%
filter(!is.na(il1_a_value))
## eGFR
model_il1a_value_adjEGFR <- lmer(log10(il1_a_value) ~ log10(egfr) +
group_sep +
(1|id), data = model_data, REML = TRUE)
Anova(model_il1a_value_adjEGFR, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il1_a_value)
## F Df Df.res Pr(>F)
## log10(egfr) 0.3496 1 447.83 0.5547
## group_sep 14.1460 9 424.92 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_il1a_value_adjEGFR)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il1_a_value) ~ log10(egfr) + group_sep + (1 | id)
## Data: model_data
##
## REML criterion at convergence: 810.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5716 -0.3349 0.0666 0.5774 3.5006
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.05197 0.228
## Residual 0.23910 0.489
## Number of obs: 500, groups: id, 186
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.59923 0.12586 4.761
## log10(egfr) 0.05827 0.09805 0.594
## group_sephealthy.000 -1.02366 0.15615 -6.556
## group_sepno_rejection.000 0.01424 0.10245 0.139
## group_sepacute_rejection.007 -0.11283 0.11920 -0.947
## group_sepno_rejection.007 -0.12177 0.11330 -1.075
## group_sepacute_rejection.090 -0.06818 0.13322 -0.512
## group_sepno_rejection.090 0.14640 0.14589 1.004
## group_sepacute_rejection.365 -0.16452 0.13886 -1.185
## group_sepno_rejection.365 0.08185 0.14728 0.556
## group_sepacute_rejection.biopsy -0.08932 0.13175 -0.678
##
## Correlation of Fixed Effects:
## (Intr) lg10() g_.000 g__.000 grp_spc_.007 grp_spn_.007
## log10(egfr) 0.730
## grp_sph.000 -0.865 -0.669
## grp_sp_.000 -0.569 0.007 0.458
## grp_spc_.007 -0.716 -0.415 0.611 0.505
## grp_spn_.007 -0.831 -0.427 0.704 0.682 0.636
## grp_spc_.090 -0.783 -0.563 0.676 0.452 0.625 0.653
## grp_spn_.090 -0.811 -0.559 0.699 0.528 0.589 0.720
## grp_spc_.365 -0.765 -0.559 0.662 0.435 0.608 0.635
## grp_spn_.365 -0.816 -0.571 0.704 0.523 0.590 0.720
## grp_spct_r. -0.637 -0.359 0.542 0.458 0.543 0.569
## grp_spc_.090 grp_spn_.090 grp_spc_.365 grp_spn_.365
## log10(egfr)
## grp_sph.000
## grp_sp_.000
## grp_spc_.007
## grp_spn_.007
## grp_spc_.090
## grp_spn_.090 0.635
## grp_spc_.365 0.652 0.621
## grp_spn_.365 0.639 0.715 0.624
## grp_spct_r. 0.556 0.524 0.537 0.525
## creatinine
model_il1a_value_adjCRE <- lmer(log10(il1_a_value) ~ log10(creatinine) +
group_sep +
(1|id), data = model_data, REML = TRUE)
Anova(model_il1a_value_adjCRE, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il1_a_value)
## F Df Df.res Pr(>F)
## log10(creatinine) 0.2582 1 438.05 0.6116
## group_sep 14.1083 9 424.95 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_il1a_value_adjCRE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il1_a_value) ~ log10(creatinine) + group_sep + (1 | id)
## Data: model_data
##
## REML criterion at convergence: 810
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5825 -0.3292 0.0701 0.5676 3.5070
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.05217 0.2284
## Residual 0.23901 0.4889
## Number of obs: 500, groups: id, 186
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.71098 0.33694 2.110
## log10(creatinine) -0.05901 0.11556 -0.511
## group_sephealthy.000 -1.01460 0.15575 -6.514
## group_sepno_rejection.000 0.01341 0.10247 0.131
## group_sepacute_rejection.007 -0.10813 0.11873 -0.911
## group_sepno_rejection.007 -0.11783 0.11341 -1.039
## group_sepacute_rejection.090 -0.06124 0.13250 -0.462
## group_sepno_rejection.090 0.15325 0.14585 1.051
## group_sepacute_rejection.365 -0.15761 0.13837 -1.139
## group_sepno_rejection.365 0.08876 0.14741 0.602
## group_sepacute_rejection.biopsy -0.08519 0.13158 -0.647
##
## Correlation of Fixed Effects:
## (Intr) lg10() g_.000 g__.000 grp_spc_.007 grp_spn_.007
## lg10(crtnn) -0.967
## grp_sph.000 -0.786 0.667
## grp_sp_.000 -0.222 0.008 0.469
## grp_spc_.007 -0.549 0.407 0.607 0.513
## grp_spn_.007 -0.608 0.429 0.705 0.688 0.635
## grp_spc_.090 -0.678 0.557 0.673 0.463 0.621 0.653
## grp_spn_.090 -0.691 0.559 0.698 0.537 0.586 0.721
## grp_spc_.365 -0.670 0.554 0.659 0.444 0.604 0.635
## grp_spn_.365 -0.702 0.572 0.704 0.531 0.587 0.721
## grp_spct_r. -0.484 0.356 0.540 0.463 0.541 0.569
## grp_spc_.090 grp_spn_.090 grp_spc_.365 grp_spn_.365
## lg10(crtnn)
## grp_sph.000
## grp_sp_.000
## grp_spc_.007
## grp_spn_.007
## grp_spc_.090
## grp_spn_.090 0.633
## grp_spc_.365 0.649 0.619
## grp_spn_.365 0.637 0.715 0.623
## grp_spct_r. 0.555 0.522 0.536 0.524
## eGFR, simpler model (biopsy and healthy people not included)
model_data <- data_long %>%
filter(!is.na(il1_a_value)) %>%
filter(group != 'healthy',
time_point != 'biopsy')
model_il1a_value_adjEGFR_simple <- lmer(log10(il1_a_value) ~ log10(egfr) +
group*time_point+
(1|id), data = model_data, REML = TRUE)
Anova(model_il1a_value_adjEGFR_simple , type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il1_a_value)
## F Df Df.res Pr(>F)
## log10(egfr) 0.5022 1 384.49 0.47897
## group 1.7232 1 124.74 0.19169
## time_point 2.4330 3 342.57 0.06483 .
## group:time_point 1.5602 3 306.31 0.19906
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_il1a_value_adjEGFR_simple)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il1_a_value) ~ log10(egfr) + group * time_point + (1 | id)
## Data: model_data
##
## REML criterion at convergence: 692.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5624 -0.2597 0.1086 0.5747 2.0997
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.05608 0.2368
## Residual 0.24028 0.4902
## Number of obs: 424, groups: id, 138
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.62772 0.11230 5.590
## log10(egfr) 0.07336 0.10298 0.712
## groupacute_rejection -0.01168 0.10338 -0.113
## time_point007 -0.14357 0.08816 -1.628
## time_point090 0.11905 0.12944 0.920
## time_point365 0.05412 0.13118 0.413
## groupacute_rejection:time_point007 0.01947 0.13016 0.150
## groupacute_rejection:time_point090 -0.20129 0.14678 -1.371
## groupacute_rejection:time_point365 -0.23553 0.15074 -1.562
##
## Correlation of Fixed Effects:
## (Intr) lg10() grpct_ tm_007 tm_090 tm_365 g_:_00 g_:_09
## log10(egfr) 0.866
## grpct_rjctn -0.279 -0.008
## time_pnt007 -0.765 -0.585 0.285
## time_pnt090 -0.754 -0.667 0.197 0.614
## time_pnt365 -0.761 -0.678 0.194 0.618 0.637
## grpct_:_007 0.172 -0.003 -0.654 -0.444 -0.150 -0.148
## grpct_:_090 0.201 0.053 -0.582 -0.229 -0.524 -0.199 0.461
## grpct_:_365 0.193 0.049 -0.568 -0.221 -0.193 -0.503 0.449 0.428There is sign that renal function affect also IL-1b, but it seems to be driven mainly by the difference of healthy subject. Analysis taking into account only planned measurements from patients does not show any indication that renal function affects the IL-1b levels.
4.5.2 IL-1b
Open code
model_data <- data_long %>%
filter(!is.na(il1_b_value))
## eGFR
model_il1a_value_adjEGFR <- lmer(log10(il1_b_value) ~ log2(egfr) +
group_sep +
(1|id), data = model_data, REML = TRUE)
Anova(model_il1a_value_adjEGFR, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il1_b_value)
## F Df Df.res Pr(>F)
## log2(egfr) 4.287 1 456.46 0.03897 *
## group_sep 61.156 9 423.86 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_il1a_value_adjEGFR)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il1_b_value) ~ log2(egfr) + group_sep + (1 | id)
## Data: model_data
##
## REML criterion at convergence: -505.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.2586 -0.4814 -0.0821 0.2832 7.0701
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.003987 0.06314
## Residual 0.015821 0.12578
## Number of obs: 500, groups: id, 186
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.21000 0.03282 36.868
## log2(egfr) 0.01602 0.00770 2.081
## group_sephealthy.000 -0.65544 0.04072 -16.094
## group_sepno_rejection.000 0.04783 0.02671 1.791
## group_sepacute_rejection.007 -0.07725 0.03074 -2.513
## group_sepno_rejection.007 -0.03984 0.02954 -1.349
## group_sepacute_rejection.090 -0.08160 0.03444 -2.369
## group_sepno_rejection.090 -0.11988 0.03798 -3.157
## group_sepacute_rejection.365 -0.10136 0.03589 -2.825
## group_sepno_rejection.365 -0.13717 0.03834 -3.577
## group_sepacute_rejection.biopsy -0.07194 0.03400 -2.116
##
## Correlation of Fixed Effects:
## (Intr) lg2(g) g_.000 g__.000 grp_spc_.007 grp_spn_.007
## log2(egfr) 0.730
## grp_sph.000 -0.865 -0.669
## grp_sp_.000 -0.568 0.007 0.457
## grp_spc_.007 -0.714 -0.420 0.609 0.497
## grp_spn_.007 -0.831 -0.427 0.704 0.688 0.632
## grp_spc_.090 -0.780 -0.569 0.675 0.444 0.627 0.649
## grp_spn_.090 -0.812 -0.560 0.700 0.533 0.587 0.725
## grp_spc_.365 -0.762 -0.563 0.660 0.428 0.610 0.631
## grp_spn_.365 -0.817 -0.571 0.704 0.528 0.588 0.725
## grp_spct_r. -0.634 -0.364 0.541 0.451 0.544 0.565
## grp_spc_.090 grp_spn_.090 grp_spc_.365 grp_spn_.365
## log2(egfr)
## grp_sph.000
## grp_sp_.000
## grp_spc_.007
## grp_spn_.007
## grp_spc_.090
## grp_spn_.090 0.634
## grp_spc_.365 0.655 0.619
## grp_spn_.365 0.637 0.721 0.622
## grp_spct_r. 0.559 0.522 0.539 0.523
## creatinine
model_il1a_value_adjCRE <- lmer(log10(il1_b_value) ~ log10(creatinine) +
group_sep +
(1|id), data = model_data, REML = TRUE)
Anova(model_il1a_value_adjCRE, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il1_b_value)
## F Df Df.res Pr(>F)
## log10(creatinine) 2.9475 1 446.93 0.0867 .
## group_sep 60.3664 9 423.90 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_il1a_value_adjCRE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il1_b_value) ~ log10(creatinine) + group_sep + (1 | id)
## Data: model_data
##
## REML criterion at convergence: -506.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.2579 -0.4712 -0.0827 0.2668 7.0904
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.004023 0.06343
## Residual 0.015848 0.12589
## Number of obs: 500, groups: id, 186
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.30700 0.08805 14.843
## log10(creatinine) -0.05210 0.03020 -1.725
## group_sephealthy.000 -0.64554 0.04068 -15.868
## group_sepno_rejection.000 0.04709 0.02675 1.760
## group_sepacute_rejection.007 -0.07221 0.03066 -2.355
## group_sepno_rejection.007 -0.03548 0.02961 -1.198
## group_sepacute_rejection.090 -0.07411 0.03429 -2.161
## group_sepno_rejection.090 -0.11234 0.03802 -2.955
## group_sepacute_rejection.365 -0.09383 0.03580 -2.621
## group_sepno_rejection.365 -0.12953 0.03843 -3.371
## group_sepacute_rejection.biopsy -0.06735 0.03399 -1.981
##
## Correlation of Fixed Effects:
## (Intr) lg10() g_.000 g__.000 grp_spc_.007 grp_spn_.007
## lg10(crtnn) -0.967
## grp_sph.000 -0.786 0.667
## grp_sp_.000 -0.222 0.008 0.468
## grp_spc_.007 -0.551 0.412 0.605 0.505
## grp_spn_.007 -0.608 0.429 0.705 0.694 0.630
## grp_spc_.090 -0.681 0.563 0.672 0.455 0.624 0.648
## grp_spn_.090 -0.691 0.559 0.699 0.542 0.584 0.726
## grp_spc_.365 -0.672 0.559 0.658 0.437 0.607 0.631
## grp_spn_.365 -0.702 0.572 0.704 0.536 0.585 0.726
## grp_spct_r. -0.487 0.361 0.539 0.456 0.542 0.564
## grp_spc_.090 grp_spn_.090 grp_spc_.365 grp_spn_.365
## lg10(crtnn)
## grp_sph.000
## grp_sp_.000
## grp_spc_.007
## grp_spn_.007
## grp_spc_.090
## grp_spn_.090 0.632
## grp_spc_.365 0.652 0.617
## grp_spn_.365 0.635 0.722 0.621
## grp_spct_r. 0.557 0.521 0.538 0.5224.5.3 IL-1RA
Open code
model_data <- data_long %>%
filter(!is.na(il1_ra_value))
## eGFR
model_il1a_value_adjEGFR <- lmer(log10(il1_ra_value) ~ log10(egfr) +
group_sep +
(1|id), data = model_data, REML = TRUE)
Anova(model_il1a_value_adjEGFR, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il1_ra_value)
## F Df Df.res Pr(>F)
## log10(egfr) 0.3599 1 470.56 0.5489
## group_sep 8.3635 9 421.80 1.647e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_il1a_value_adjEGFR)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il1_ra_value) ~ log10(egfr) + group_sep + (1 | id)
## Data: model_data
##
## REML criterion at convergence: 328.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.3925 -0.6091 -0.1129 0.4015 3.4088
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.02745 0.1657
## Residual 0.08422 0.2902
## Number of obs: 500, groups: id, 186
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.14221 0.07774 40.419
## log10(egfr) -0.03650 0.06055 -0.603
## group_sephealthy.000 -0.54051 0.09652 -5.600
## group_sepno_rejection.000 0.11687 0.06332 1.846
## group_sepacute_rejection.007 -0.10337 0.07129 -1.450
## group_sepno_rejection.007 -0.04941 0.07003 -0.706
## group_sepacute_rejection.090 -0.08403 0.08022 -1.047
## group_sepno_rejection.090 -0.03606 0.08969 -0.402
## group_sepacute_rejection.365 -0.25929 0.08351 -3.105
## group_sepno_rejection.365 -0.14709 0.09055 -1.624
## group_sepacute_rejection.biopsy -0.15415 0.07893 -1.953
##
## Correlation of Fixed Effects:
## (Intr) lg10() g_.000 g__.000 grp_spc_.007 grp_spn_.007
## log10(egfr) 0.730
## grp_sph.000 -0.864 -0.669
## grp_sp_.000 -0.568 0.007 0.457
## grp_spc_.007 -0.708 -0.429 0.605 0.482
## grp_spn_.007 -0.830 -0.427 0.703 0.699 0.622
## grp_spc_.090 -0.776 -0.579 0.671 0.429 0.631 0.639
## grp_spn_.090 -0.814 -0.560 0.700 0.544 0.583 0.735
## grp_spc_.365 -0.757 -0.571 0.656 0.414 0.614 0.622
## grp_spn_.365 -0.818 -0.571 0.705 0.539 0.584 0.735
## grp_spct_r. -0.630 -0.373 0.538 0.437 0.548 0.556
## grp_spc_.090 grp_spn_.090 grp_spc_.365 grp_spn_.365
## log10(egfr)
## grp_sph.000
## grp_sp_.000
## grp_spc_.007
## grp_spn_.007
## grp_spc_.090
## grp_spn_.090 0.630
## grp_spc_.365 0.660 0.615
## grp_spn_.365 0.634 0.734 0.619
## grp_spct_r. 0.563 0.519 0.543 0.520
## creatinine
model_il1a_value_adjCRE <- lmer(log10(il1_ra_value) ~ log10(creatinine) +
group_sep +
(1|id), data = model_data, REML = TRUE)
Anova(model_il1a_value_adjCRE, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il1_ra_value)
## F Df Df.res Pr(>F)
## log10(creatinine) 0.0000 1 462.62 0.9966
## group_sep 9.0465 9 421.81 1.543e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_il1a_value_adjCRE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il1_ra_value) ~ log10(creatinine) + group_sep + (1 | id)
## Data: model_data
##
## REML criterion at convergence: 328.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.3829 -0.5912 -0.1104 0.4121 3.4478
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.02792 0.1671
## Residual 0.08401 0.2898
## Number of obs: 500, groups: id, 186
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.1772369 0.2085549 15.235
## log10(creatinine) -0.0003088 0.0715359 -0.004
## group_sephealthy.000 -0.5796413 0.0964005 -6.013
## group_sepno_rejection.000 0.1171901 0.0633890 1.849
## group_sepacute_rejection.007 -0.1218307 0.0709551 -1.717
## group_sepno_rejection.007 -0.0675076 0.0701644 -0.962
## group_sepacute_rejection.090 -0.1121288 0.0797432 -1.406
## group_sepno_rejection.090 -0.0665761 0.0897262 -0.742
## group_sepacute_rejection.365 -0.2881525 0.0831712 -3.465
## group_sepno_rejection.365 -0.1785399 0.0906953 -1.969
## group_sepacute_rejection.biopsy -0.1720801 0.0787602 -2.185
##
## Correlation of Fixed Effects:
## (Intr) lg10() g_.000 g__.000 grp_spc_.007 grp_spn_.007
## lg10(crtnn) -0.967
## grp_sph.000 -0.786 0.667
## grp_sp_.000 -0.221 0.008 0.468
## grp_spc_.007 -0.556 0.422 0.601 0.489
## grp_spn_.007 -0.608 0.429 0.704 0.706 0.620
## grp_spc_.090 -0.686 0.573 0.668 0.439 0.628 0.638
## grp_spn_.090 -0.692 0.560 0.700 0.554 0.579 0.736
## grp_spc_.365 -0.676 0.568 0.654 0.422 0.611 0.621
## grp_spn_.365 -0.703 0.572 0.705 0.548 0.581 0.736
## grp_spct_r. -0.491 0.370 0.535 0.442 0.546 0.555
## grp_spc_.090 grp_spn_.090 grp_spc_.365 grp_spn_.365
## lg10(crtnn)
## grp_sph.000
## grp_sp_.000
## grp_spc_.007
## grp_spn_.007
## grp_spc_.090
## grp_spn_.090 0.628
## grp_spc_.365 0.657 0.613
## grp_spn_.365 0.632 0.735 0.617
## grp_spct_r. 0.561 0.517 0.541 0.5194.5.4 IL-36b
Open code
model_data <- data_long %>%
filter(!is.na(il36_b_value))
## eGFR
model_il1a_value_adjEGFR <- lmer(log10(il36_b_value) ~ log10(egfr) +
group_sep +
(1|id), data = model_data, REML = TRUE)
Anova(model_il1a_value_adjEGFR, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il36_b_value)
## F Df Df.res Pr(>F)
## log10(egfr) 1.3490 1 483.16 0.246
## group_sep 8.4876 9 419.12 1.086e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_il1a_value_adjEGFR)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il36_b_value) ~ log10(egfr) + group_sep + (1 | id)
## Data: model_data
##
## REML criterion at convergence: 568.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.3808 -0.3381 0.0098 0.4113 3.1354
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.05718 0.2391
## Residual 0.13051 0.3613
## Number of obs: 500, groups: id, 186
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.79421 0.10017 7.929
## log10(egfr) 0.09076 0.07778 1.167
## group_sephealthy.000 -0.61739 0.12452 -4.958
## group_sepno_rejection.000 0.05236 0.08192 0.639
## group_sepacute_rejection.007 -0.21499 0.08931 -2.407
## group_sepno_rejection.007 -0.09652 0.09049 -1.067
## group_sepacute_rejection.090 -0.43762 0.10105 -4.331
## group_sepno_rejection.090 -0.21276 0.11517 -1.847
## group_sepacute_rejection.365 -0.60102 0.10508 -5.720
## group_sepno_rejection.365 -0.26771 0.11628 -2.302
## group_sepacute_rejection.biopsy -0.55169 0.09899 -5.573
##
## Correlation of Fixed Effects:
## (Intr) lg10() g_.000 g__.000 grp_spc_.007 grp_spn_.007
## log10(egfr) 0.728
## grp_sph.000 -0.863 -0.666
## grp_sp_.000 -0.570 0.006 0.458
## grp_spc_.007 -0.700 -0.439 0.599 0.462
## grp_spn_.007 -0.829 -0.425 0.701 0.715 0.608
## grp_spc_.090 -0.768 -0.591 0.665 0.409 0.637 0.625
## grp_spn_.090 -0.815 -0.559 0.701 0.561 0.576 0.748
## grp_spc_.365 -0.750 -0.581 0.650 0.395 0.620 0.608
## grp_spn_.365 -0.820 -0.570 0.705 0.555 0.578 0.748
## grp_spct_r. -0.624 -0.384 0.532 0.418 0.552 0.544
## grp_spc_.090 grp_spn_.090 grp_spc_.365 grp_spn_.365
## log10(egfr)
## grp_sph.000
## grp_sp_.000
## grp_spc_.007
## grp_spn_.007
## grp_spc_.090
## grp_spn_.090 0.624
## grp_spc_.365 0.666 0.609
## grp_spn_.365 0.628 0.750 0.613
## grp_spct_r. 0.568 0.514 0.548 0.515
## creatinine
model_il1a_value_adjCRE <- lmer(log10(il36_b_value) ~ log10(creatinine) +
group_sep +
(1|id), data = model_data, REML = TRUE)
Anova(model_il1a_value_adjCRE, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
##
## Response: log10(il36_b_value)
## F Df Df.res Pr(>F)
## log10(creatinine) 1.4245 1 476.83 0.2333
## group_sep 8.5998 9 419.39 7.347e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_il1a_value_adjCRE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il36_b_value) ~ log10(creatinine) + group_sep + (1 | id)
## Data: model_data
##
## REML criterion at convergence: 567.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.3756 -0.3412 0.0117 0.4202 3.1371
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.05725 0.2393
## Residual 0.13044 0.3612
## Number of obs: 500, groups: id, 186
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.01992 0.26817 3.803
## log10(creatinine) -0.11026 0.09196 -1.199
## group_sephealthy.000 -0.61969 0.12431 -4.985
## group_sepno_rejection.000 0.05102 0.08192 0.623
## group_sepacute_rejection.007 -0.21531 0.08895 -2.421
## group_sepno_rejection.007 -0.09802 0.09058 -1.082
## group_sepacute_rejection.090 -0.43842 0.10049 -4.363
## group_sepno_rejection.090 -0.21481 0.11516 -1.865
## group_sepacute_rejection.365 -0.60224 0.10471 -5.752
## group_sepno_rejection.365 -0.27011 0.11640 -2.320
## group_sepacute_rejection.biopsy -0.55244 0.09882 -5.590
##
## Correlation of Fixed Effects:
## (Intr) lg10() g_.000 g__.000 grp_spc_.007 grp_spn_.007
## lg10(crtnn) -0.967
## grp_sph.000 -0.784 0.665
## grp_sp_.000 -0.222 0.008 0.468
## grp_spc_.007 -0.561 0.432 0.595 0.470
## grp_spn_.007 -0.607 0.427 0.702 0.720 0.606
## grp_spc_.090 -0.692 0.585 0.663 0.420 0.634 0.625
## grp_spn_.090 -0.693 0.559 0.701 0.568 0.573 0.749
## grp_spc_.365 -0.681 0.578 0.648 0.404 0.616 0.608
## grp_spn_.365 -0.703 0.571 0.706 0.562 0.575 0.749
## grp_spct_r. -0.497 0.381 0.531 0.424 0.550 0.543
## grp_spc_.090 grp_spn_.090 grp_spc_.365 grp_spn_.365
## lg10(crtnn)
## grp_sph.000
## grp_sp_.000
## grp_spc_.007
## grp_spn_.007
## grp_spc_.090
## grp_spn_.090 0.622
## grp_spc_.365 0.664 0.607
## grp_spn_.365 0.627 0.750 0.612
## grp_spct_r. 0.566 0.512 0.546 0.5145 IL18 forms plot
At first, we need to transform IL-18 levels to molar concentration (pM/L)
Open code
data_long <- data_long %>%
mutate( il18bp_pM = (((il18_bp_value/1e12)*1000) / 19000)*1e12,
il18tot_pM = (((il18_value/1e12)*1000) / 17200)*1e12,
il18free_pM = (((il18_free/1e12)*1000) / 17200) *1e12)Now, we can create plots to demonstrate the relationship between IL-18BP levels and IL-18 levels (both free and total). We will concentrate on healthy individuals and the initial measurement from the group without rejection, noted for its significant variability in IL-18BP and IL-18 levels. The plots will display examples showcasing both minimal and very high IL-18 binding.
Open code
## prepare data
il18d <- data_long %>%
filter(!is.na(il18_value),
!is.na(il18_bp_value),
!is.na(il18_free)) %>%
mutate(group = factor(group, levels =
c("no_rejection",
"acute_rejection",
"healthy")),
time_point = factor(time_point)) %>%
mutate(time_point = if_else(group == 'healthy', 'healthy', time_point)) %>%
mutate(group_sep = interaction(group, time_point))
## Plot for healthy
p18_healthy <- il18d %>%
filter(group_sep == 'healthy.healthy') %>%
ggplot(aes(x = il18bp_pM)) +
geom_point(aes(y = il18free_pM, color = "IL-18 Free")) +
geom_point(aes(y = il18tot_pM, color = "IL-18 Total"), pch=1) +
geom_line(data = . %>%
select(il18bp_pM, il18free_pM, il18tot_pM) %>%
pivot_longer(cols = -il18bp_pM,
names_to = "variable", values_to = "value"),
aes(x = il18bp_pM, y = value, group = il18bp_pM), col = 'grey20') +
scale_color_manual(values = c("IL-18 Free" = "blue", "IL-18 Total" = "red")) +
labs(x = "IL-18BP (pM)", y = "IL-18 (pM)", color = "") +
coord_cartesian(ylim = c(0, NA)) +
ggtitle("Samples from healthy donors")
p18_healthy
Open code
if(file.exists('gitignore/figures/p18_healthy.pdf') == FALSE){
ggsave("gitignore/figures/p18_healthy.pdf",
plot = p18_healthy, width = 8, height = 5.5, units = "in")
}
## plot for D0, group without rejection
p18_no_reject_0 <- il18d %>%
filter(group_sep == 'no_rejection.000') %>%
ggplot(aes(x = il18bp_pM)) +
geom_point(aes(y = il18free_pM, color = "IL-18 Free")) +
geom_point(aes(y = il18tot_pM, color = "IL-18 Total"), pch=1) +
geom_line(data = . %>%
select(il18bp_pM, il18free_pM, il18tot_pM) %>%
pivot_longer(cols = -il18bp_pM,
names_to = "variable", values_to = "value"),
aes(x = il18bp_pM, y = value, group = il18bp_pM), col = 'grey20') +
scale_color_manual(values = c("IL-18 Free" = "blue", "IL-18 Total" = "red")) +
labs(x = "IL-18BP (pM)", y = "IL-18 (pM)", color = "") +
coord_cartesian(ylim = c(0, NA)) +
ggtitle("No rejection group, the 1st measurement")
p18_no_reject_0
Open code
if(file.exists('gitignore/figures/p18_no_reject_0.pdf') == FALSE){
ggsave("gitignore/figures/p18_no_reject_0.pdf",
plot = p18_no_reject_0, width = 8, height = 5.5, units = "in")
}
## plot for D0, group with rejection
p18_acute_reject_0 <- il18d %>%
filter(group_sep == 'acute_rejection.000') %>%
ggplot(aes(x = il18bp_pM)) +
geom_point(aes(y = il18free_pM, color = "IL-18 Free")) +
geom_point(aes(y = il18tot_pM, color = "IL-18 Total"), pch=1) +
geom_line(data = . %>%
select(il18bp_pM, il18free_pM, il18tot_pM) %>%
pivot_longer(cols = -il18bp_pM,
names_to = "variable", values_to = "value"),
aes(x = il18bp_pM, y = value, group = il18bp_pM), col = 'grey20') +
scale_color_manual(values = c("IL-18 Free" = "blue", "IL-18 Total" = "red")) +
labs(x = "IL-18BP (pM)", y = "IL-18 (pM)", color = "") +
coord_cartesian(ylim = c(0, NA)) +
ggtitle("Acute rejection group, the 1st measurement")
p18_acute_reject_0
Open code
if(file.exists('gitignore/figures/p18_acute_reject_0.pdf') == FALSE){
ggsave("gitignore/figures/p18_acute_reject_0.pdf",
plot = p18_acute_reject_0, width = 8, height = 5.5, units = "in")
}6 Reproducibility
Open code
sessionInfo()
## R version 4.3.2 (2023-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=cs_CZ.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=cs_CZ.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=cs_CZ.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=cs_CZ.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/Prague
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] lqmm_1.5.8 effects_4.2-2 rstudioapi_0.13 quantreg_5.88
## [5] SparseM_1.81 coxed_0.3.3 survival_3.5-7 rms_6.7-1
## [9] Hmisc_5.1-1 bayesplot_1.8.1 ggdist_3.3.0 kableExtra_1.3.4
## [13] lubridate_1.8.0 corrplot_0.92 arm_1.13-1 lme4_1.1-34
## [17] MASS_7.3-60 janitor_2.2.0 RJDBC_0.2-10 rJava_1.0-6
## [21] DBI_1.1.2 projpred_2.6.0 brms_2.16.3 Rcpp_1.0.11
## [25] glmnet_4.1-8 Matrix_1.6-3 boot_1.3-28 cowplot_1.1.1
## [29] pROC_1.18.0 mgcv_1.9-1 nlme_3.1-163 openxlsx_4.2.5
## [33] flextable_0.9.3 sjPlot_2.8.15 car_3.1-2 carData_3.0-5
## [37] gtsummary_1.7.2 emmeans_1.7.2 ggpubr_0.4.0 forcats_1.0.0
## [41] stringr_1.5.0 dplyr_1.1.3 purrr_1.0.2 readr_2.1.2
## [45] tidyr_1.2.0 tibble_3.2.1 ggplot2_3.4.3 tidyverse_1.3.1
##
## loaded via a namespace (and not attached):
## [1] fs_1.6.3 matrixStats_1.0.0 httr_1.4.2
## [4] webshot_0.5.5 insight_0.19.5 SparseGrid_0.8.2
## [7] tools_4.3.2 backports_1.4.1 sjlabelled_1.2.0
## [10] utf8_1.2.3 R6_2.5.1 DT_0.17
## [13] withr_2.5.0 Brobdingnag_1.2-7 prettyunits_1.1.1
## [16] gridExtra_2.3 cli_3.6.1 textshaping_0.3.6
## [19] performance_0.10.5 gt_0.9.0 shinyjs_2.1.0
## [22] officer_0.6.2 sandwich_3.0-1 labeling_0.4.2
## [25] sass_0.4.7 mvtnorm_1.1-3 polspline_1.1.24
## [28] ggridges_0.5.3 askpass_1.1 commonmark_1.9.0
## [31] systemfonts_1.0.4 StanHeaders_2.21.0-7 foreign_0.8-86
## [34] gfonts_0.2.0 svglite_2.1.0 readxl_1.3.1
## [37] httpcode_0.3.0 generics_0.1.2 shape_1.4.6
## [40] gtools_3.9.2 crosstalk_1.2.0 distributional_0.3.2
## [43] zip_2.2.0 inline_0.3.19 loo_2.4.1
## [46] fansi_1.0.4 abind_1.4-5 lifecycle_1.0.3
## [49] multcomp_1.4-18 yaml_2.3.5 snakecase_0.11.1
## [52] grid_4.3.2 promises_1.2.0.1 crayon_1.5.2
## [55] miniUI_0.1.1.1 lattice_0.22-5 haven_2.4.3
## [58] pillar_1.9.0 knitr_1.37 estimability_1.3
## [61] shinystan_2.5.0 codetools_0.2-19 glue_1.6.2
## [64] fontLiberation_0.1.0 data.table_1.14.8 broom.helpers_1.14.0
## [67] vctrs_0.6.3 cellranger_1.1.0 gtable_0.3.4
## [70] assertthat_0.2.1 xfun_0.40 mime_0.12
## [73] survey_4.2-1 coda_0.19-4 rsconnect_0.8.25
## [76] iterators_1.0.14 shinythemes_1.2.0 ellipsis_0.3.2
## [79] TH.data_1.1-0 pbkrtest_0.5.2 xts_0.12.1
## [82] fontquiver_0.2.1 threejs_0.3.3 rstan_2.21.3
## [85] tensorA_0.36.2 rpart_4.1.23 colorspace_2.1-0
## [88] nnet_7.3-19 tidyselect_1.2.0 processx_3.8.2
## [91] compiler_4.3.2 curl_4.3.2 rvest_1.0.2
## [94] htmlTable_2.4.0 xml2_1.3.3 fontBitstreamVera_0.1.1
## [97] bayestestR_0.13.1 colourpicker_1.1.1 posterior_1.2.0
## [100] checkmate_2.0.0 scales_1.2.1 dygraphs_1.1.1.6
## [103] callr_3.7.3 digest_0.6.29 minqa_1.2.4
## [106] rmarkdown_2.25 htmltools_0.5.6 pkgconfig_2.0.3
## [109] base64enc_0.1-3 highr_0.9 dbplyr_2.1.1
## [112] fastmap_1.1.0 rlang_1.1.1 htmlwidgets_1.6.2
## [115] shiny_1.5.0 farver_2.1.1 zoo_1.8-9
## [118] jsonlite_1.7.3 magrittr_2.0.3 Formula_1.2-4
## [121] munsell_0.5.0 gdtools_0.3.3 stringi_1.7.12
## [124] plyr_1.8.8 pkgbuild_1.3.1 parallel_4.3.2
## [127] sjmisc_2.8.9 ggeffects_1.3.1 splines_4.3.2
## [130] hms_1.1.1 sjstats_0.18.2 ps_1.6.0
## [133] igraph_1.5.1 uuid_1.0-3 ggsignif_0.6.3
## [136] markdown_1.8 reshape2_1.4.4 stats4_4.3.2
## [139] rstantools_2.1.1 crul_1.4.0 reprex_2.0.1
## [142] evaluate_0.15 mitools_2.4 RcppParallel_5.1.7
## [145] modelr_0.1.8 nloptr_2.0.0 tzdb_0.4.0
## [148] foreach_1.5.2 httpuv_1.6.5 MatrixModels_0.5-0
## [151] openssl_1.4.6 broom_1.0.5 xtable_1.8-4
## [154] rstatix_0.7.0 later_1.3.0 viridisLite_0.4.2
## [157] ragg_1.2.1 cluster_2.1.6 gamm4_0.2-6
## [160] bridgesampling_1.1-2